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1.  Accomplishments 

1.1.  Summary  of  accomplishments 

The  mission  of  the  Defense  Threat  Reduction  Agency  requires  the  quantitative  study  and  accu¬ 
rate  prediction  for  complex  multiphysics  systems  that  couple  together  physical  processes  spanning 
wide  range  of  scales  in  behavior.  Treatment  of  such  systems  depends  on  accurate  numerical  sim¬ 
ulation  of  mathematical  models  expressed  as  systems  of  partial  differential  equations  posed  on 
domains  with  complicated  geometry.  Prediction  of  the  behavior  involves  treating  the  propagation 
of  stochastic  uncertainty  through  the  mathematical  models  and  solving  inverse  problems  for  de¬ 
termining  parameters  based  on  observations  on  model  output.  Quantifying  the  accuracy  of  such 
computations  requires  accurate  estimation  of  the  numerical  error  in  quantities  of  interest  com¬ 
puted  from  numerical  solutions  that  take  into  account  all  sources  of  error,  e.g.  from  discretization, 
representation  of  geometry,  finite  sampling. 

This  project  focuses  on  development  of  mathematical  tools  for  dealing  with  these  problems  in 
the  context  of  multiphysics  models  of  interest  using  relevant  numerical  methods  to  the  mission  of 
the  DTRA.  The  main  approach  is  a  posteriori  error  analysis  based  on  computable  residuals,  solu¬ 
tion  of  adjoint  problems,  and  variational  analysis.  This  approach  estimates  the  error  in  specified 
quantities  of  interest.  Computable  residuals  involving  the  approximate  solution  are  used  to  quan¬ 
tify  the  size  of  various  discretization  errors  while  the  solution  of  adjoint  equations  (generalized 
Green’s  functions)  are  used  to  quantify  the  effects  of  stability  in  producing  errors.  Much  of  the 
project  dealt  with  dealing  the  significant  mathematical  issues  that  arise  when  numerically  solv¬ 
ing  complex  multiphysics  models.  Practical  computational  constraints  requires  the  use  of  a  wide 
variety  of  discretization  approaches,  e.g.  operator  decomposition  and  splitting,  explicit  time  inte¬ 
gration,  iterative  solution  methods  with  few  iterations,  finite  volume  and  specialized  finite  differ¬ 
ence  methods.  The  introduction  of  such  techniques  complicates  both  the  identification  of  suitable 
residuals  and  definition  of  suitable  adjoint  problems',  The  project  also  dealt  with  issues  arising  in 
“multi-discretization”  approaches,  when  various  components  of  a  coupled  system  are  solved  with 
different  numerical  methods  and  numerical  grids.  Another  focus  was  the  treatment  of  problems 
posed  on  complex  domains,  e.g.  on  manifold  surfaces  in  space  and/or  on  domains  with  complex 
boundaries.  In  this  case,  the  goal  was  to  treat  the  effects  of  inaccuracies  and/or  uncertainty  in  the 
representation  of  the  domain  geometry.  Finally,  we  also  establshed  several  rigorous  convergence 
results  for  a  class  of  goal-oriented  adaptive  methods  that  are  designed  to  drivdriving  the  error  in  a 
specific  quantity  of  interest  below  a  given  tolerance. 

Along  with  theoretical  development,  the  project  studied  the  practical  implementation  of  a  pos¬ 
teriori  error  estimates  for  complex  physics,  including  high  performance  issues.  The  project  also 
addressed  the  question  of  efficient  computation.  The  availability  of  accurate  error  estimates  raises 
the  ability  to  develop  efficient  adaptive  error  control  algorithms  in  which  various  discretization 
parameters  are  adjusted  based  on  relative  contributions  to  the  overall  error  in  order  to  achieve  a 
desired  accuracy  with  minimal  computational  work.  In  another  direct,  the  project  expanded  a  pos¬ 
teriori  error  estimates  for  computed  distributions  and  probabilities  arising  in  computational  sensi¬ 
tivity  analysis  and  developed  generalized  adaptive  algorithms  that  allow  for  balancing  all  sources 
of  error  and  uncertainty  affecting  the  analysis. 

The  project  P.I.s’  undertook  a  significant  degree  of  interdisciplinary  interaction  during  the 
projects  in  order  to  insure  that  project  accomplishments  would  have  impact  in  science  and  en¬ 
gineering. 

1 .2.  Detailed  descriptions  of  specific  accomplishments 

In  this  section,  we  describe  specific  technical  accomplishments  of  the  project. 

A  posteriori  error  analysis  for  a  transient  conjugate  heat  transfer 

We  analyzed  the  accuracy  of  an  operator  decomposition  finite  element  method  for  a  transient  con¬ 
jugate  heat  transfer  problem  consisting  of  two  materials  coupled  through  a  common  boundary.  We 
derive  accurate  a  posteriori  error  estimates  that  account  for  the  transfer  of  error  between  compo- 


3 


nenls  of  the  operator  decomposition  method  as  well  as  the  errors  in  solving  the  iterative  system. 
We  address  a  loss  of  order  of  convergence  that  results  from  the  decomposition,  and  show  that  the 
order  of  convergence  is  limited  by  the  accuracy  of  the  transferred  gradient  information.  We  ex¬ 
tend  a  boundary  flux  recovery  method  to  transient  problems  and  use  it  to  regain  the  expected  order 
of  accuracy  in  an  efficient  manner.  In  addition,  we  use  the  a  posteriori  error  estimates  to  adap¬ 
tively  compute  the  recovered  boundary  flux  only  within  the  domain  of  dependence  for  a  quantity 
of  interest. 

A  posteriori  error  estimation  and  adaptive  mesh  refinement  for  a  multiscale  operator  decomposi¬ 
tion  approach  to  fluid-solid  heat  transfer 

We  analyze  a  multiscale  operator  decomposition  finite  element  method  for  a  conjugate  heat  transfer 
problem  consisting  of  a  fluid  and  a  solid  coupled  through  a  common  boundary.  We  derive  accurate 
a  posteriori  error  estimates  that  account  for  all  sources  of  error,  and  in  particular  the  transfer  of  error 
between  fluid  and  solid  domains.  We  use  these  estimates  to  guide  adaptive  mesh  refinement.  In 
addition,  we  provide  compelling  numerical  evidence  that  the  order  of  convergence  of  the  operator 
decomposition  method  is  limited  by  the  accuracy  of  the  transferred  gradient  information,  and  adapt 
a  so-called  boundary  flux  recovery  method  developed  for  elliptic  problems  in  order  to  regain  the 
optimal  order  of  accuracy  in  an  efficient  manner.  In  an  appendix,  we  provide  an  argument  that 
explains  the  numerical  results  provided  sufficient  smoothness  is  assumed. 

Nonparametric  density  estimation  for  randomly  perturbed  elliptic  problems 

We  study  the  nonparametric  density  estimation  problem  for  a  quantity  of  interest  computed  from 
solutions  of  an  elliptic  partial  differential  equation  with  randomly  perturbed  coefficients  and  data. 
We  derive  an  efficient  method  for  computing  samples  and  generating  an  approximate  probability 
distribution  based  on  Lion’s  domain  decomposition  method  and  the  Neumann  series.  We  then 
derive  an  a  posteriori  error  estimate  for  the  computed  probability  distribution  reflecting  all  sources 
of  deterministic  and  statistical  errors.  Finally,  we  develop  an  adaptive  error  control  algorithm 
based  on  the  a  posteriori  estimate,  we  extend  the  analysis  to  include  a  “modeling  error”  term  that 
accounts  for  the  effects  of  the  resolution  of  the  statistical  description  of  the  random  variation  and 
modify  the  adaptive  algorithm  to  adapt  the  resolution  of  the  statistical  description.  We  also  prove 
some  related  convergence  results. 

A  posteriori  error  analysis  for  cell-centered  finite  volume  methods  for  semilinear  elliptic  problems 

We  conduct  a  goal-oriented  a  posteriori  analysis  for  the  error  in  a  quantity  of  interest  computed 
from  a  cell-centered  finite  volume  scheme  for  a  semilinear  elliptic  problem.  To  carry  out  the 
analysis,  we  use  an  equivalence  between  the  cell-centered  finite  volume  scheme  and  a  mixed  finite 
element  method  with  special  choice  of  quadrature. 

Blockwise  adaptivity  for  time  dependent  problems  based  on  coarse  scale  adjoint  solutions 

We  describe  and  test  an  adaptive  algorithm  for  evolution  problems  that  employs  a  sequence  of 
“blocks”  consisting  of  fixed,  though  non-uniform,  space  meshes.  This  approach  offers  the  advan¬ 
tages  of  adaptive  mesh  refinement  but  with  reduced  overhead  costs  associated  with  load  balancing, 
re-meshing,  matrix  reassembly,  and  the  solution  of, adjoint  problems  used  to  estimate  discretiza¬ 
tion  error  and  the  effects  of  mesh  changes.  We  describe  several  strategies  to  determine  appropriate 
block  discretizations  from  coarse  scale  solution  information  using  adjoint-based  a  posteriori  error 
estimates  and  demonstrate  the  behavior  of  the  algorithms  in  a  set  of  examples. 

Conservative  discretization  and  a  posteriori  error  analysis  for  a  cut  cell  diffusion  problems  with 
complex  geometry 

We  study  the  solution  of  a  diffusive  process  in  a  domain  where  the  diffusion  coefficient  changes 
discontinuously  across  a  curved  interface.  We  consider  discretizations  that  use  regularly-shaped 
meshes,  so  that  the  interface  “cuts”  through  the  cells  (elements  or  volumes)  without  respecting 
the  regular  geometry  of  the  mesh.  Consequently,  the  discontinuity  in  the  diffusion  coefficients  has 
a  strong  impact  on  the  accuracy  and  convergence  of  the  numerical  method.  This  motivates  the 
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derivation  of  computational  error  estimates  that  yield  accurate  estimates  for  specified  quantities  of 
interest.  For  this  purpose,  we  adapt  the  well-known  adjoint  based  a  posteriori  error  analysis  tech¬ 
nique  used  for  finite  element  methods.  In  order  to  employ  this  method,  we  describe  a  systematic 
approach  to  discretizing  a  cut-cell  problem  that  handles  complex  geometry  in  the  interface  in  a 
natural  fashion  yet  reduces  to  the  well-known  Ghost  Fluid  Method  in  simple  cases.  We  test  the 
accuracy  of  the  estimates  in  a  series  of  examples. 

A  measure -theoretic  computational  method  for  inverse  sensitivity  problems 

We  consider  the  inverse  sensitivity  analysis  problem  of  quantifying  the  uncertainty  of  inputs  to  a 
deterministic  map  given  specified  uncertainty  in  a  linear  functional  of  the  output  of  the  map.  This 
is  a  version  of  the  model  calibration  or  parameter  estimation  problem  for  a  deterministic  map.  We 
assume  that  the  uncertainty  in  the  quantity  of  interest  is  represented  by  a  random  variable  with 
a  given  distribution  and  we  use  the  Law  of  Total  Probability  to  express  the  inverse  problem  for 
the  corresponding  probability  measure  on  the  input  space.  Assuming  that  the  map  from  the  input 
space  to  the  quantity  of  interest  is  smooth,  we  solve  the  generally  ill-posed  inverse  problem  by 
using  the  Implicit  Function  Theorem  to  derive  a  method  for  approximating  the  set-valued  inverse 
that  provides  an  approximate  quotient  space  representation  of  the  input  space.  We  then  derive  an 
efficient  computational  approach  to  compute  a  measure  theoretic  approximation  of  the  probability 
measure  on  the  input  space  imparted  by  the  approximate  set-valued  inverse  that  solves  the  inverse 
problem.  We  also  treat  the  situation  in  which  the  output  of  the  map  is  determined  implicitly  and 
is  difficult  and/or  expensive  to  evaluate,  e.g  requiring  the  solution  of  a  differential  equation,  and 
hence  the  output  of  the  map  is  approximated  numerically.  The  main  goal  is  an  a  posteriori  error 
estimate  that  can  be  used  to  evaluate  the  accuracy  of  the  computed  distribution  solving  the  inverse 
problem  taking  into  account  all  sources  of  statistical  and  numerical  deterministic  errors.  We  present 
a  general  analysis  for  the  method  and  then  apply  the  analysis  to  the  case  of  a  map  determined  by 
the  solution  of  an  initial  value  problem. 

A  posteriori  analysis  of  multirate  numerical  methods  for  multiscale  ordinary  differential  equations 

We  analyze  a  multirate  time  integration  method  for  systems  of  ordinary  differential  equations  that 
present  significantly  different  scales  within  the  components  of  the  model.  We  interpret  the  mul¬ 
tirate  method  as  a  multiscale  operator  decomposition  method  and  use  this  formulation  to  conduct 
both  an  a  priori  error  analysis  and  a  hybrid  a  priori  ~  a  posteriori  error  analysis.  The  hybrid  analy¬ 
sis  has  the  form  of  a  computable  a  posteriori  leading  order  expression  and  a  provably-higher  order 
a  priori  expression.  Both  analyses  distinguish  the  effects  of  the  discretization  of  each  component 
from  the  effects  of  multirate  solution.  The  effects  on  stability  arising  from  the  multirate  solution 
are  reflected  in  perturbations  to  certain  associated  adjoint  operators. 

Convergence  theory  for  goal-oriented  adaptive  methods 

In  the  first  of  the  convergence  theory  subprojects  of  the  DTRA  project.  We  developed  a  new  con¬ 
vergence  theory  for  a  general  class  of  adaptive  approximation  algorithms  for  nonlinear  operator 
equations,  and  then  used  the  theory  to  obtain  convergence,  contraction,  and  optimality  results  for 
practical  adaptive  finite  element  methods  (AFEM)  applied  to  several  classes  of  nonlinear  elliptic 
equations  and  systems  of  elliptic  equations.  The  results  can  be  viewed  as  extending  the  recent  con¬ 
vergence  results  for  linear  problems  of  Morin,  Siebert  and  Veeser,  and  of  Nochetto  et.  al  to  more 
general  nonlinear  problems  (with  G.  Tsogtgerel  and  Y.  Zhu).  We  also  develop  new  mathematical 
results  for  hierarchical  error  indicators  to  drive  AFEM  algorithms,  and  establish  condition  number 
estimates  for  appropriate  preconditioners  (with  J.  Ovall  and  R.  Szypowski).  We  have  further  ex¬ 
tended  these  results  to  the  class  of  adaptive  methods  that  were  the  target  of  this  DTRA  research 
probejet:  goal-oriented  adaptive  methods  that  are  designed  to  drive  the  error  in  a  quantity  of  inter¬ 
est  below  a  given  tolerance.  In  2009,  Mommer  and  Stevenson  developed  a  goal-oriented  adaptive 
method  for  the  Poisson  equation,  together  with  rigorous  convergence  and  complexity  results  for 
their  method,  establishing  what  was  apparently  the  first  convergence  result  for  a  goal-oriented 
adaptive  method.  We  have  now  extended  the  results  of  Mommer  and  Stevenson  to  goal-oriented 
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adaptive  methods  for  general  linear  convection-diffusion  elliptic  problems  (with  S.  Pollock).  In  a 
second  manuscript,  these  results  were  further  extended  to  a  large  class  of  scalar  nonlinear  prob¬ 
lems  (with  S.  Pollock  and  Y.  Zhu).  All  three  articles  have  now  been  posted  on  arXiv,  submitted 
for  publication,  and  are  currently  in  review.  All  of  the  techniques  are  demonstrated  for  practical 
problems  of  interest  using  the  FETK  software  (see  below). 

Analysis  of  multiphysics  problems  with  complex  domains 

We  analyzed  a  large  class  of  regularized  Navier-Stokes  and  Magnetohydrodynamics  (MHD)  mod¬ 
els  in  three-dimensional  spatial  domains,  a  class  which  includes  the  Navier-Stokes  equations,  the 
Navier-Stokes-alpha  model,  the  Leray-alpha  model,  the  Modified  Leray-alpha  model,  the  Simpli¬ 
fied  Bardina  model,  the  Navier-Stokes- Voight  model,  the  Navier-Stokes-alpha-like  models,  and 
certain  MHD  models,  in  addition  to  representing  a^  larger'3-parameter  family  of  models  not  previ¬ 
ously  analyzed.  We  recovered  a  number  of  known  results  for  established  models,  but  also  obtained 
new  results  for  all  models  in  this  general  family,  including  existence,  regularity,  uniqueness,  sta¬ 
bility,  attractor  existence  and  dimension,  and  existence  of  determining  operators.  (J.  Nonlinear 
Science  2009,  with  E.  Lunasin  and  G.  Tsogtgerel.) 

We  then  develop  and  analyze  numerical  methods  for  approximation  of  stationary  and  evolution 
problems  on  surfaces,  including  coupled  elliptic-parabolic  systems.  A  major  theoretical  break- 
trough  was  showing  how  the  recent  finite  element  error  estimates  of  Demlow  and  Dziuk  can  be 
recovered  from  a  more  general  approach  involving  the  analysis  of  variational  crimes  in  Hilbert 
complexes,  generalizing  their  results  for  surface  finite  elements  to  arbitrary  spatial  dimension  and 
to  applications  involving  higher-dimensional  differential  forms  and  both  linear  and  nonlinear  equa¬ 
tions.  This  generalization  was  made  possible  through  the  use  and  extension  of  finite  element  exte¬ 
rior  calculus  (FEEC).  (Found.  Comput.  Math.  2012,  with  A.  Stem.)  We  have  now  extended  this 
work  in  FEEC  in  the  direction  of  time-dependent  problems;  we  completed  and  submitted  a  new 
manuscript  in  2012  that  extends  these  results  on  surface  finite  element  methods  to  scalar  parabolic 
and  hyperbolic  problems,  including  again  nonlinear  problems  (with  A.  Gillette).  We  also  give  an 
analysis  of  the  singularities  in  a  fundamentally  important  model  in  biochemistry,  and  develop  a 
number  of  AFEM-based  numerical  techniques  for  treating  these  degenerate  features  in  a  provably 
high-fidelity  way  (Comm.  Comput.  Phys.  2012,  with  J.  McCammon,  Y.  Zhou,  Y.  Zhu,  Z.  Yu). 

In  addition,  we  have  developed  and  implemented  goal-oriented,  adjoint-based,  a  posteriori 
error  estimates  for  elliptic  problems  on  smooth  manifolds.  In  particular,  the  estimates  take  into 
account  the  effects  of  domain  curvature  on  accuracy.  We  also  considered  the  problem  of  small 
random  perturbations  to  the  manifold,  pointing  the  way  to  treat  problems  in  which  the  domain 
is  determined  experimentally  or  by  measurement.  This  work  is  nearing  completion  and  will  be 
submitted  in  Summer  2012  (with  W.  Newton) 

Analysis  of  elliptic  problems  on  domains  with  randomly  perturbed  boundaries 

We  developed  a  systematic  approach  to  solve  elliptic  problems  on  domains  that  have  randomly 
perturbed  boundaries,  after  first  classifying  such  problems  into  .several  different  classes.  The  results 
are  particularly  relevant  to  situations  in  which  the  boundaries  are  obtained  through  measurement 
or  are  subject  to  error.  The  approach  avoids  the  need  to  remesh  each  new  domain  in  a  random 
sampling  Monte  Carlo  solution.  Moreover,  we  derive  a  posteriori  error  estimates  that  indicate  how 
random  perturbations  in  the  boundary  affect  the  accuracy  of  computed  solutions. 

A  posteriori  error  analysis  of  explicit,  IMEX,  and  truncated  Picard  iteration  time  integration  meth¬ 
ods 

Explicit,  Implicit/Explicit  (IMEX),  and  truncated  Picard  iteration  time  integration  methods  are 
widely  employed  to  solve  multiphysics  applications  in  defense  and  department  of  energy  enter¬ 
prises,  e.g.  such  as  reacting  flows.  Such  methods  requires  significant  alterations  for  a  posteriori 
error  analysis  in  order  to  describe  the  effects  of  these  approaches  on  both  stability  and  accuracy. 
Therefore,  last  year  we  undertook  the  systematic  study  of  a  posteriori  error  analysis  for  explicit, 
truncated  Picard  iteration,  and  implicit/explicit  (IMEX)  time  integration  methods.  For  explicit 
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methods,  we  introduce  special  projection  operators  into  the  standard  finite  element  formulation  for 
evolution  problems.  These  projection  operators  are  (1)  a  truncated  Taylor  expansion  computed  at 
a  past  time  node  and  (2)  extrapolation  from  a  interpolatory  polynomial  using  values  at  a  collection 
of  previous  nodes.  We  then  alter  the  a  posteriori  error  analysis  to  include  terms  that  measure  the  ef¬ 
fects  of  these  projections,  yielding  distinct  “explicit”  time  integration  terms  in  the  a  posteriori  error 
analysis.  We  recently  have  extended  this  approach  to  treat  IMEX  methods.  To  analyze  truncated 
Picard  iteration  methods,  we  exploit  an  old  result  of  H.  Keller  and  J.  Keller  for  the  “matricant”, 
which  is  the  exponential  form  of  the  solution  operator  of  a  linear  non-autonomous  evolution  prob¬ 
lem.  This  provides  a  way  to  define  the  adjoint  for  a  solution  obtained  by  truncated  Picard  iteration, 
which  we  then  use  in  the  a  posteriori  error  analysis.  We  have  also  extended  this  analysis  to  implicit 
methods  that  employ  Jacobi  iteration  to  solve  the  systems  at  each  step. 

Coupled  parabolic-elliptic  systems 

Estep  and  Holst  collaborated  on  the  development  methods  and  a  posteriori  error  analysis  for  cou¬ 
pled  parabolic-elliptic  systems  of  equations.  The  main  application  is  on  modeling  of  black  holes. 
A  new  development  in  the  Holst  group  has  been  the  extension  of  their  recent  work  on  finite  ele¬ 
ment  exterior  calculus  to  parabolic  and  hyperbolic  problems  (completed  and  submitted  in  2012), 
which  will  provide  a  very  strong  mathematical  framework  for  the  development  of  methods  and  a 
posteriori  analysis  for  coupled  parabolic-elliptic  problems.  This  extension  to  FEEC  is  now  being 
combined  with  our  recent  work  on  goal-oriented  adaptive  methods  using  a  variational  framework, 
by  which  the  elliptic  component  of  the  system  is  combined  with  implicit  time-stepping  schemes 
to  provide  constraints  in  a  Lagrange  multiplier  formulation.  We  are  able  to  show  convergence 
for  the  adaptive  scheme,  generalizing  our  recent  work  on  convergence  theory  for  goal-oriented 
adaptive  methods  (with  S.  Pollock,  Y.  Zhu). 

Coupled  ordinary  differential  equation  -  parabolic  differential  equation 

Estep  and  Hameed  (along  with  collaborators)  derived  and  implemented  a  posteriori  error  estimates 
for  systems  of  evolution  equations  consisting  of ’a  reaction-diffusion  problem  posed  on  a  global 
domain  coupled  to  systems  of  ordinary  differential  equations  in  a  collection  of  small  cells  parti¬ 
tioning  the  global  domain.  The  local  cell  problems  model  chemical  reactions  that  determine  the 
local  physical  conditions  driving  the  parabolic  problem.  The  analysis  takes  into  account  the  itera¬ 
tion  error  in  solving  the  coupled  systems. 

New  approaches  to  adaptive  error  control  for  evolution  problems 

Estep  and  Hameed  (along  with  collaborators)  developed  new  adaptive  error  control  algorithms  that 
take  into  account  cancellation  of  errors  to  improve  efficiency.  The  approach  identifies  periods  of 
time  over  which  there  is  significant  cancellation.  Inside  the  regions,  uniform  refinements  are  used 
to  preserve  the  favorable  cancellation,  while  the  time  step  sizes  in  the  various  regions  are  adjusted 
according  to  the  contribution  to  the  overall  error  from  the  regions. 

Implementation  of  theoretical  results 

The  last  major  goal  in  this  project  is  implementation  of  the  theoretical  results  into  the  FETK  code. 
For  this  purpose,  we  recruited  a  full  time  postdoc,  Ryan  Szypowski,  working  at  UCSD  under  the 
supervision  of  co-PI  Michael  Holst  with  responsibility  to  carry  out  the  implementation  and  testing. 
He  is  being  jointly  supervised  by  the  PI  D.  Estep.  This  FETK  deveopment  has  focused  on  providing 
a  robust,  theory-based  convergent  adaptive  finite  element  implementation  for  nonlinear  problems 
which  retains  linear  complexity.  This  has  included  work  on  the  following  specific  components, 
which  have  been  implemented  in  both  the  MATLAB  subset  FETKLab  of  the  2D  code  in  FETK  as 
well  as  in  the  full  2D/3D  code  in  FETK: 

)  . 

1.  The  clement  marking  strategy  was  updated  to  be  based  on  “Dorller  Marking”.  Special  care 
was  taken  to  use  a  linear-time  complexity  binning  approach  as  opposed  to  an  actual  sort. 
Only  this  type  of  marking  strategy,  which  is  not  often  used  in  practice  due  to  its  poten¬ 
tial  costs  unless  carefully  implemented,  allows  for  establishing  both  convergence  and  linear 
overall  computational  complexity  of  the  adaptive  algorithms. 
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2.  A  number  of  new  error  estimators  were  added.  They  include: 

(a)  A  hierarchical  error  estimator  based  on  face-bump  functions  which  was  proven  in  our 
recent  publications  to  be  efficient,  reliable,  and  robust.  This  work  included  the  addition 
of  a  new  cubic  bump  finite  element  space,  which  led  to  a  better  understanding  of  how 
we  can  improve  the  finite  element  space  implementation  to  allow  for  future  additions. 

(b)  An  error  estimator  based  on  the  solution  of  a  dual  problem,  which  we  refer  to  as  dual- 
weighted  residual  (DWR).  This  implementation  involved  leveraging  the  work  on  the 
bump-function  library  above,  as  well  as  the  development  of  high-order  quadrature 
rules,  and  the  ability  to  maintain  two  distinct  and  unrelated  adaptive  meshes  during 
a  computation,  with  quantities  being  projected  back  and  forth  between  the  meshes  as 
needed. 

(c)  An  error  estimator  based  on  smoothed  gradients.  This  is  based  on  recent  work  of  R. 
Bank  and  J.  Xu,  collaborators  of  the  Pis. 

3.  W.  Newton,  co-advised  by  Estep  and  Holst,  implemented  the  a  posteriori  error  estimates  that 
account  for  error  in  the  description  of  the  manifold  on  which  the  problem  is  posed  developed 
in  his  thesis. 

4.  A  driver  application  for  solving  nonlinear  problems  using  inexact  Newton  solvers  based  on 
a  multilevel  approach  was  written.  This  has  been  used  for  most  of  the  problems  described 
above. 

5.  Prior  to  2013,  FETK  and  the  FETKLab  MATLAB  subset  of  FETK  were  primarily  based  on 
linear  finite  element  discretizations,  with  enough  partial  support  for  higher-order  elements 
to  allow  for  the  use  of  e.g.  bump  functions  in  error  indicators  and  formulation  of  dual  prob¬ 
lems.  A  general  element  class  was  developed  in  early  2013  to  allow  for  use  of  any  type 
of  Lagrange-type  element  for  either  the  primarl  or  dual  problem.  Both  linear  and  quadratic 
elements  were  then  implemented  and  are  provided  with  the  FETK  code  base  as  element 
examples.  Our  recent  manuscripts  with  new  convergence  results  for  goal-oriented  methods 
contain  a  large  collection  of  numerical  examples  that  now  exploit  this  infrastructure  to  care¬ 
fully  compare  a  number  of  adaptive  methods  based  on  goal  functions  (with  S.  Pollock  and 
Y.  Zhu). 

2.  Training  and  Professional  Development 

The  support  of  this  project  has  partially  cohtributed  to  the  training  and  professional  devel¬ 
opment  for  three  graduate  students  and  three  postdocs.  This  includes  specialized  research-level 
instruction  and  individual  mentoring  as  well  as  participation  in  large  research  group  activities  di¬ 
rected  by  the  Pis.  Students  and  postdocs  were  encouraged  to  participate  in  professional  meetings 
and  to  interact  with  researchers  in  other  universities  and  in  national  DOE  laboratories  as  appropri¬ 
ate.  Students  and  postdocs  were  trained  to  write  and  prepare  and  deliver  professional  presentations. 

Details  for  the  trainees: 

•  Will  Newton  received  his  Ph.D.  from  CSU  in  2011,  and  then  was  hired  as  a  Research  Scien¬ 
tist  Class  1  in  PI  Estep’s  group.  His  primary  focus  is  a  project  on  multiscale  models  of  new 
nuclear  fuels  supported  by  a  contract  from  Idaho  National  Laboratory.  He  has  continued  to 
work  on  research  related  to  this  project  following  up  on  the  work  in  his  thesis.  Thesis  is  “A 
Posteriori  Error  Estimates  for  the  Poisson  Problem  on  Closed,  Two-Dimensional  Surfaces”, 
available  from  Colorado  State  University  Library. 

•  Nate  Burch  received  his  Ph.D.  from  CSU  in  2011,  and  then  took  a  two  year  postdoc  po¬ 
sition  at  SAMSI  (Statistical  and  Mathematical  Sciences  Institute)  as  part  of  the  Program 
on  Uncertainty  Quantification.  Thesis  is  “Probabilistic  Foundation  of  Nonlocal  Diffusion 
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and  Formulation  and  Analysis  for  Elliptic  Problems  on  Uncertain  Domains”,  available  from 
Colorado  State  University  Library. 

•  The  CSU  postdoc  Jehanzeb  Hamecd  is  in  the  second  year  of  his  position  in  PI  Estep’s  group. 
His  primary  focus  is  a  project  on  a  Department  of  Energy  Uncertainty  Quantification  project 
that  is  jointly  conducted  with  Sandia  National  Laboratory.  Part  of  his  research  is  related  to 
the  activities  supported  in  this  project. 

•  Jonny  Serenesa  received  his  Ph.D.  from  UCSD  in  2012,  and  has  been  doing  pre-  and  post¬ 
doctoral  work  at  UC  Davis.  His  doctoral  work  was  jointly  supervised  by  PI  Holst  and  S. 
Shkoller  at  UC  Davis,  and  he  is  currently  working  for  a  startup  company  in  the  Bay  Area. 

•  Ryan  Szypowski  received  his  Ph.D.  from  UCSD  in  2008,  and  remained  at  UCSD  working 
with  Holst  as  a  postdoc  and  then  research  scientist  until  2012.  He  moved  to  a  tenure-track 
position  in  the  Mathematics  Department  at  Cal  Poly  Pomona  in  Fall  2012. 

•  Andrew  Gillette  received  his  Ph.D.  from  UT  Austin  in  2011,  and  joined  Holst’s  group  at 
UCSD  as  a  postdoctoral  fellow  in  Fall  2011.  He  helped  push  forward  both  the  the  project 
involving  Ryan  Szypowski,  and  the  development  of  an  FEEC-based  error  analysis  frame¬ 
work  for  parabolic  and  hyperbolic  problems.  In  Fall  2013,  Andrew  is  starting  a  tenure-track 
faculty  position  in  the  mathematics  department  at  the  University  of  Arizona. 

•  Sara  Pollock  received  her  Ph.D.  from  UCSD  2012,  and  remained  at  UCSD  working  with 
Holst  as  a  postdoc  during  the  2012-2012  academic  year.  In  Fall  2013,  Sara  is  starting  a 
3-year  named  postdoctoral  position  in  the  mathematics  department  at  Texas  A&M. 

3.  Dissemination  ' 

We  have  disseminated  the  research  in  this  project  through  submission  of  peer-reviewed  research 
articles,  presenfing  many  invited  talks  at  universities  and  conferences,  and  publishing  software 
developed  in  this  project  for  public  access.  A  summary  of  this  activity  during  this  project: 

•  53  research  articles  related  to  the  project  research  have  appeared  or  are  accepted 

•  1 9  research  articles  related  to  the  project  research  are  currently  under  review 

•  5  book  and/or  book  chapters  have  appeared  or  are  being  written 

•  60  invited  lectures  at  universities  and  professional  meetings 

Applications  to  multiscale/multiphysics  physical  and  engineering  systems 

In  conjunction  with  collaborators  in  engineering,  chemistry  and  biophysics,  we  have  applied  many 
of  the  algorithms  and  techniques  for  multiphysics  and  multiscale  problems  developed  in  this 
DTRA-supported  research  program.  Our  focus  continues  to  be  on  applications  in  material,  chemi¬ 
cal  and  biological  physics  of  relevance  to  DOD,  DTRA,  and  DOE  missions.  In  addition  to  our  pub¬ 
lications  placed  in  the  mathematics  literature,  we  have  placed  joint  publications  from  these  research 
collaborations  with  physical  scientists  and  engineers  in  a  broad  spectrum  of  leading  scientific  jour¬ 
nals  to  maximize  the  impact  of  our  results,  including:  Physical  Review  Letters,  Physical  Review 
D,  Journal  of  Nonlinear  Science,  Classical  and  Quantum  Gravity,  Journal  of  Chemical  Theory 
and  Computation,  Journal  of  Cell  Science,  Journal  of  Structural  Biology,  Biophysical  Journal, 
PLoS  Computational  Biology,  IMA  Journal  on  Applied  Mathematics,  Computer  Aided  Geomet¬ 
ric  Design,  BIT,  Applied  Numerical  Mathematics,  IEEE  Journal  on  Engineering  in  Medicine  and 
Biology,  IEEE  Transactions  on  Biomedical  Computing,  Frontiers  in  Computational  Physiology 
and  Medicine,  Investigative  Ophthalmology  and  Visual  Science,  Journal  of  Scientific  Computing, 
Journal  of  Applied  Mathematics  and  Computation,  Communications  in  Computational  Physics, 
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Journal  of  Molecular  Graphics  and  Modeling,  Journal  of  Physical  Chemistry  B,  Journal  of  Chem¬ 
ical  Physics,  Communications  in  Mathematical  Physics,  Annals  of  Nuclear  Engineering,  Journal 
of  Computational  Physics,  Acta  Biomaterialia,  Computer  Methods  in  Applied  Mechanics  and  En¬ 
gineering,  Journal  of  Engineering  Mathematics,  and  Foundations  of  Computational  Mathematics. 

4.  Products 

4.1.  Publications,  conference  papers,  and  presentations 

The  following  papers  were  accepted  or  appeared  during  March  27,  2009  -  September  1,  2009 

•  A  posteriori  analysis  and  adaptive  error  control  for  multiscale  operator  decomposition  meth¬ 
ods  for  coupled  elliptic  systems  I:  One  way  coupled  systems,  V.  Carey,  D.  Estep,  and  S. 
Tavener,  SIAM  Journal  on  Numerical  Analysis  47  (2009),  740-761 

•  A  posteriori  error  analysis  for  a  transient  conjugate  heat  transfer  problem,  D.  Estep,  S. 
Tavener,  T.  Wildey,  Finite  Elements  in  Analysis  and  Design  45  (2009),  263-271 

•  Nonparametric  density  estimation  for  randomly  perturbed  elliptic  problems  I:  Computa¬ 
tional  methods,  a  posteriori  analysis,  and  adaptive  error  control,  D.  Estep,  A.  Malqvist,  and 
S.  Tavener,  SIAM  Journal  on  Scientific  Computing  31  (2009),  2935-2959 

•  Solving  the  Einstein  constraints  on  multi-block  triangulations  using  finite  elements,  O.  Ko- 
robkin,  B.  Aksoylu,  M.  Holst,  E.  Pazos,  and  M.  Tiglio,  Class.  Quant.  Grav.  26  (2009),  No. 
14,  145007  (28  pp).  (arXiv:gr-qc/0801.1823) 

•  An  adaptive  finite  element  method  for  solving  the  exact  Kohn-Sham  equation  of  density  func¬ 
tional  theory,  E.  Bylaska,  M.  Holst,  and  J.  Weare,  Journal  of  Chemical  Theory  and  Compu¬ 
tation,  5  (2009),  pp.  937-948. 

•  Finite  Element  Analysis  of  Drug  Electrostatic  Diffusion:  Inhibition  Rate  Studies  in  N1  Neu¬ 
raminidase,  Y.  Cheng,  M.  Holst,  and  J.A.  McCammon,  Biocomputing  2009:  Proceedings  of 
the  Pacific  Symposium,  R.B.  Altman,  A.K.  Dunker,  L.  Hunter,  T.  Murray,  and  T.E.  Klein, 
eds.,  2009,  pp.  281-292. 

,  ■  'I 

•  Three-dimensional  reconstruction  reveals  new  details  of  membrane  systems  for  calcium  sig¬ 
naling  in  the  heart,  T.  Hayashi,  M.E.  Martone,  Z.  Yu,  A.  Thor,  M.  Doi,  M.  Holst,  M.H. 
Elli.sman,  and  M.  Hoshijima,  J.  Cell  Sci.,  Vol.  122  (April,  2009),  No.  7,  pp.  1005-1013. 

•  Rough  Solutions  of  the  Einstein  Constraints  on  closed  manifolds  without  near-CMC  condi¬ 
tions,  M.  Holst,  G.  Nagy,  and  G.  Tsogtgerel,  Comm.  Math.  Phys.,  Vol.  288  (June  2009), 
No.  2,  pp.  547-613.  (arXiv;gr-qc/07 12.0798) 

•  Multi-Scale  Modeling  of  Ventricular  Myocytes:  Contributions  of  structural  and  functional 
heterogeneities  to  excitation-contraction  coupling  in  the  normal  and  failing  rodent  heart,  S. 
Lu,  A.  Michailova,  J.  Saucerman,  Y.  Cheng  Z.  Yu,  T.  Kaiser,  W.  Li,  R.  Bank,  M.  Holst,  A. 
McCammon,  T.  Hayashi,  M.  Hoshijima,  P.  Arzbcrger,  and  A.  McCulloch,  IEEE  Journal  on 
Engineering  in  Medicine  and  Biology,  Vol.  28  (March-April  2009),  No.  2,  pp.  46-57. 

•  Convergence  and  Optimality  of  Adaptive  Mixed  Finite  Element  Methods,  L.  Chen,  M.  Holst, 
and  J.  Xu,  Math.  Comp.,  Vol.  78  (2009),  No.  265,  pp.  33-53. 
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The  following  papers  were  accepted  or  appeared  during  September  2,  2009  -  September  1,  2010 

•  Nonparametric  density  estimation  for  randomly  perturbed  elliptic  problems  II:  Applications 
and  adaptive  modeling,  D.  Estep,  A.  Malqvist,  S.  Tavener,  International  Journal  for  Numer¬ 
ical  Methods  in  Engineering  80  (2009),  846-867 

•  A  posteriori  error  analysis  of  a  cell-centered  finite  volume  method  for  semilinear  elliptic 
problems,  D.  Estep,  M.  Pemice,  D.  Pham,  S.  Tavener,  H.  Wang,  Journal  of  Computational 
and  Applied  Mathematics  233  (2009),  459  -  472 

•  A  posteriori  error  estimation  and  adaptive  mesh  refinement  for  a  multi-discretization  oper¬ 
ator  decomposition  approach  to  fluid-solid  heat  transfer,  D.  Estep,  S.  Tavener,  T  Wildey, 
Journal  of  Computational  Physics  229  (2010),  4143-4158 

•  Blockwise  adaptivity  for  time  dependent  problems  based  on  coarse  scale  adjoint  solutions, 
V.  Carey,  D.  Estep,  A.  Johansson,  M.  Larson,  and  S.  Tavener,  SIAM  Journal  on  Scientific 
Computing  32  (2010),  2121  -  2145 

•  Numerical  analysis  of  Ca2-\-  signaling  in  rat  ventricular  myocytes  with  realistic  transverse- 
axial  tubular  geometry  and  inhibited  sarcoplasmic  reticulum,  Y.  Cheng,  Z.  Yu,  M.  Hoshi- 
jima,  M.  Holst,  A.  McCulloch,  and  J.  M.  ad  A.P.  Michailova,  PLoS  Computational  Biology 
6(2010),  pp.  el 000972: 1-16. 

•  Poisson-Nernst-Planck  equations  for  simulation  biomolecular  diffusion-reaction  processes 
I:  Finite  element  solutions,  B.  Lu,  M.  Holst,  J.  McCammon,  and  Y.  Zhou,  J.  of  Comput 
Phys.  229  (2010),  6679-7794  (16  pp). 

•  Analysis  of  a  general  family  of  regularized  Navier-Stokes  and  MHD  models,  M.  Holst,  E. 
Lunasin,  and  G.  Tsogtgerel,  J.  Nonlin.  Sci.,  20  (2010),  pp.  523-567. 

The  following  book  chapter  appeared  during  September  2,  2009  -  September  1,  2010 

•  Error  estimation  for  multiscale  operator  decomposition  for  multiphysics  problems,  D.  Es¬ 
tep,  Chapter  1 1 ,  in  Bridging  the  Scales  in  Science  and  Engineering,  J.  Fi.sh,  editor,  Oxford 
University  Press,  2010 

The  following  books  were  under  contract  or  appeared  during  September  2,  2009  -  April  5,2013 

•  Practical  Analysis  in  Many  Variables,  D.  Estep,  SIAM,  2010. 

•  Green’s  Functions  and  Boundary  Value  Problems,  Third  Edition,  I.  Stakgold  and  M.  Holst, 
John-Wiley,  888  pages,  February  2011. 

.1,1 

The  following  nonrefereed  papers  appeared  during  September  2,  2009  -  September  1,  2010 

•  CSE  2009:  Graduate  Education  in  CSE  -  Structure  for  the  Zoo?,  H.-J.  Bungartz  and  D. 
Estep,  SIAM  News  42,  2009 

•  Computational  Science  and  Engineering  Education:  SIAM’S  Perspective,  H.-J.  Bungartz,  D. 
Estep,  U.  Rude,  and  P.  Turner,  IEEE  Computing  in  Science  and  Engineering  1 1  (2009),  5-11 

•  Interview  with  Chief  Editor  of  the  SIAM  CSE  Book  Series,  D.  Estep,  SIAM  News  43  (2010) 


11 


The  following  papers  were  accepted  or  appeared  during  September  2,  2010  -  September  1,  201 1 

•  A  computational  measure  theoretic  method  for  inverse  sensitivity  problems  I:  Basic  method 
and  analysis,  J.  Breidt,  T.  Butler,  and  D.  Estep,  SIAM  Journal  on  Numerical  Analysis,  201 1, 
49  (2011),  1836-1859 

•  /I  posteriori  error  analysis  for  a  cut  cell  finite  volume  method,  D.  Estep,  S.  Tavener,  M. 
Pemice,  H.  Wang,  Computer  Methods  in  Applied  Mechanics  and  Engineering,  2010,  233 
(2009),  459-472 

•  Parameter  estimation  and  directional  leverage  with  applications  in  differential  equations,  N. 
Burch,  D.  Estep,  and  J.  Hoeting,  Metrica,  Metrika,  DOI:  10.1007/s00184-01 1-0358-4,  201 1 

•  Continuum  Modeling  and  Control  of  Large  Mobile  Networks,  Y.  Zhang,  E.  K.  P.  Chong,  J. 
Hannig,  and  D.  Estep,  Proceedings  of  the  49th  Annual  Allerton  Conference  on  Communica¬ 
tion,  Control  and  Computing,  Illinois,  201 1 

•  Nonparameteric  density  estimation  for  randomly  perturbed  elliptic  problems  111:  Conver¬ 
gence.  complexity,  and  generalizations,  D.  Estep,  M.  Holst,  and  A.  Malqvist,  Journal  of 
Applied  Mathematics  and  Computing  38  (2012),  367-387 

•  An  efficient,  reliable  and  robust  error  estimator  for  elliptic  problems  in  R^,  M.  Holst,  J. 
Ovall,  and  R.  Szypowski,  Applied  Numerical  Mathematics,  61  (201 1),  675695 

•  Efficient  mesh  optimization  schemes  based  on  optimal  delaunay  triangulations,  L.  Chen  and 
M.  Holst,  Computer  Methods  in  Applied  Mechanics  and  Engineering  200  (201 1),  967984 

•  Adaptive  finite  element  modeling  techniques  for  the  Poisson-Boltzmann  equation,  M.  Holst, 
J.  McCammon,  Z.  Yu,  Y.  Zhou,  and  Y.  Zhu,  Communications  in  Computational  Physics  1 1 
(2012),  pp.  179-214. 

•  Convergence  analysis  of  finite  element  approximations  of  the  Joule  heating  problem  in  three 
spatial  dimensions,  M.  Holst,  M.  Larson,  A.  Malqvist,  and  R.  Soderlund,  BIT  50  (2010) 
pp. 781-795. 

•  Semilinear  mixed  problems  on  Hilbert  complexes  and  their  numerical  approximation,  M. 
Holst  AND  A.  Stem,  Foundations  of  Computational  Mathematics,  2010, 12  (2012),  pp.  363- 
387 

•  Adaptive  solution  of  the  Poisson-Boltzmann  equation  using  goal-oriented  error  indicators, 
B.  Aksoylu,  S.  Bond,  E.  Cyr,  AND  M.  Holst,  J.  Sci.  Comput.  52  (2012),  202-225  (23  pp).  ’ 

The  following  papers  were  accepted  or  appeared  during  September  2,  2011  -  September  1,  2012 

•  A  computational  measure  theoretic  approach  to  inverse  sensitivity  problems  II:  A  posteriori 
error  analysis,  T.  Butler,  D.  Estep  and  J.  Sandelin,  SIAM  Journal  on  Numerical  Analysis  50 
(2012)  ^ 

•  Viscoelastic  Effects  During  Loading  Play  an  Integral  Role  in  Soft  Tissue  Mechanics,  K. 
Troyer,  D.  Estep,  and  C.  Puttlitz,  Acta  Biomaterialia  8  (2012),  234-244 

•  A  posteriori  analysis  of  multirate  numerical  method  for  ordinary  differential  equations,  D. 
Estep,  V.  Ginting,  S.  Tavener,  2012,  Computer  Methods  in  Applied  Mechanics  and  Engi¬ 
neering,  223-224  (2012),  10-27 

•  Adaptive  error  control  for  an  elliptic  optimization  problem,  App\\cdib\Q  \m\ys\s  D  Estep 
and  S.  Lee,  2012,  D01:10.1080/0003681 1.2012.683785,  1-15 
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•  Analysis  of  routing  protocols  and  interference-limited  communication  in  large  networks  via 
continuum  modeling,  N.  Burch,  E.  Chong,  D.  Estep,  J.  Hannig,  Journal  of  Engineering  Math¬ 
ematics,  2012,  (DOI)  10. 1007/s  10665-0 12-9566-9 

•  A  numerical  method  for  solving  a  stochastic  inverse  problem  for  parameters,  T.  Butler  and 
D.  Estep,  Annals  of  Nuclear  Energy,  2012,  10.1016/j.anucene.2012.05.016 

•  Geometric  variational  crimes:  Hilbert  complexes,  finite  element  exterior  calculus,  and  prob¬ 
lems  on  hypersurfaces,  M.  Holst  and  A.  Stem,  Foundations  of  Computational  Mathematics, 
12(2012),  pp.  263-293. 

•  Multi-scale  modeling  of  calcium  dynamics  in  ventricular  myocytes  with  realistic  transverse 
tubules,  Z.  Yu,  G.  Yao,  M.  Hoshijima,  A.  Michailova,  and  M.  Holst,  IEEE  TBME  Let¬ 
ters,  Special  Issue  on  Multi-Scale  Modeling  and  Analysis  for  Computational  Biology  and 
Medicine,  58  (2011),  No.  10,  2947-2951  (4  pp). 

•  Multiscale  continuum  modeling  and  simulation  of  biological  processes:  From  molecular 
electro-diffusion  to  sub-cellular  signaling  transduction,  Y.  Cheng,  M.  Holst,  J.  McCammon, 
and  A.  Michailova,  Comput.  Sci.  Disc.,  5  (2012),  015002-015015  (13  pp). 

•  The  Navier-Stokes-Voight  model  for  image  inpainting,  M.  Ebrahimi,  M.  Holst,  and  E.  Lu- 
nasin,  IMA  J.  Appl.  Math.,  doi:10.1093/imamatyhxr069  (2012),  1-26  (26  pp). 

•  Numerical  bifurcation  analysis  of  conformal  formulations  of  the  Einstein  constraints,  M.  Holst 
and  V.  Kungurtsev,  Phys.  Rev.  D,  84  (2011),  pp.  124038(1)-! 24038(8). 

•  Modeling  cardiac  calcium  sparks  in  a  three-dimensional  reconstruction  of  a  calcium  re¬ 
lease  unit,  J.  Hake,  A.  Edwards,  Z.  Yu,  P.  Kekenes-Huskey,  A.  Michailova,  A.  McCammon, 
M.  Holst,  M.  Hoshijima,  and  A.  McCulloch,  J.  Physiol.,  590  (2012),  No.  18, 4403-4422  (18 
pp)- 

•  Localized  glaucomatous  change  detection  within  the  proper  orthogonal  decomposition  frame¬ 
work,  M.  Balasubramanian,  D.  Kriegman,  C.  Bowd,  M.  Holst,  R.  Winreb,  P.  Sample,  and 
L.  Zangwill,  Invest.  Ophthalmol.  Vis.  Sci.,  53  (2012),  No.  7,  3615-3628  (14  pp). 

•  Quality  tetrahedral  mesh  smoothing  via  boundary-optimized  Delaunay  triangulation,  Z.  Gao, 
Z.  Yu,  and  M.  Holst,  Computer  Aided  Geometric  Design,  29(9):707-721,  2012. 

•  Modeling  effects  of  L-type  Ca2+  current  and  Na-\r-Ca2+  exchanger  on  Ca2+  trigger  flux  in 
rabbit  myocytes  with  realistic  T-tubule  geometries,  P.  Kekenes-Huskey,  Y.  Cheng,  J.  Hake, 

F.  Sachse,  J.  Bridge,  M.  Holst,  J.  McCammon,  A.  McCulloch,  and  A.  Michailova,  Frontiers 
in  Physiology,  3  (2012),  pp.  1-14. 

The  following  papers  were  accepted,  appeared  or  were  submitted  and  still  pending  review  during 
September  2,  2011  -  September  1,  2012 

•  A  Posteriori  Analysis  and  Adaptive  Error  Control  for  Multiscale  Operator  Decomposition 
Solution  of  Elliptic  Systems  11:  Fully  Coupled  Systems,  V.  Carey,  D.  Estep,  S.  Tavener,  Inter¬ 
national  Journal  of  Numerical  Methods  in  Engineering,  201 1,  in  revision 

•  A  posteriori  analysis  of  an  iterative  multi-discretization  method  for  reaction-diffusion  sys¬ 
tems,  J.  H.  Chaudhry,  D.  Estep,  V.  Ginting,  and  S.  Tavener,  Computer  Methods  in  Applied 
Mechanics  and  Engineering,  2012,  in  revision 

•  A-posteriori  error  estimates  for  mixed  finite  element  and  finite  volume  methods  for  problems 
coupled  through  a  boundary  with  non-matching  grids,  T.  Arbogast,  D.  Estep,  B.  Sheehan, 
and  S.  Tavener,  IMA  J.  Numerical  Analysis,  2012,  in  revision 
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•  Multilevel  preconditioners  for  discontinuous  Galerkin  approximations  of  elliptic  problems 
with  jump  coefficients,  B.  Ayuso  de  Dios,  M.  Holst,  Y.  Zhu,  and  L.  Zikatanov,  in  Proceedings 
of  the  Twentieth  International  Conference  on  Domain  Decomposition  Methods,  San  Diego, 
USA,  San  Diego,  CA,  USA,  February  201 1. 

t  .  ‘ 

•  Local  multilevel  preconditioners  for  elliptic  equations  with  jump  coefficients  on  bisection 
grids,  L.  Chen,  M.  Holst,  J.  Xu,  and  Y.  Zhu,  Submitted  for  publication. 

•  Local  convergence  of  adaptive  methods  for  nonlinear  partial  differential  equations,  M.  Holst, 

G.  Tsogtgercl,  and  Y  Zhu,  Submitted  for  publication. 

•  The  Lichnerowicz  equation  on  compact  manifolds  with  boundary,  M.  Holst  and  G.  Tsogt- 
gerel.  Submitted  for  publication. 

•  Adaptive  finite  element  methods  with  inexact  solvers  for  the  nonlinear  Poisson-Boltzmann 
equation,  M.  Holst,  R.  Szypowski,  and  Y.  Zhu,  in  Proceedings  of  the  Twentieth  International 
Conference  on  Domain  Decomposition  Methods,  San  Diego,  USA,  San  Diego,  CA,  USA, 
February  2011. 

•  Barrier  methods  for  critical  exponent  problems  in  geometric  analysis  and  mathematical 
physics,  J.  Erway  and  M.  Holst,  Submitted  for  publication. 

•  Finite  element  error  estimates  for  critical  exponent  semilinear  problems  without  angle  con¬ 
ditions,  R.  Bank,  M.  Holst,  R.  Szypowski,  and  Y.  Zhu,  Submitted  for  publication. 

•  Convergence  and  optimality  of  goal-orientied  adaptive  finite  element  methods  for  nonsym- 
metric  problems,  M.  Holst  and  S.  Pollock,  Submitted  for  publication. 

•  Generalized  solutions  to  semilinear  elliptic  PDF  with  applications  to  the  Lichnerowicz  equa¬ 
tion,  M.  Holst  and  C.  Meier,  Submitted  for  publication. 

•  Finite  element  exterior  calculus  for  evolution  problems,  A.  Gillette  and  M.  Holst,  Submitted 
for  publication. 

•  Two-grid  methods  for  semilinear  interface  problems,  M.  Holst,  R.  Szypowski,  and  Y.  Zhu, 
Accepted  for  publication  in  Numer.  Methods  Partial  Differtial  Equations. 

•  Convergence  of  goal-oriented  adaptive  finite  element  methods  for  semilinear  problems,  M.  Holst, 

S.  Pollock,  and  Y.  Zhu,  Submitted  for  publication. 

•  Feature-preserving  surface  mesh  smoothing  via  suboptional  Delaunay  triangulation,  Z.  Gao, 

Z.  Yu,  and  M.  Holst,  Graphical  Models,  75  (2013),  pp.  23-38. 

The  following  papers  were  accepted,  appeared  or  were  submitted  and  still  pending  review  during 
September  2,  2012  -  April  5,  2012 

•  Multiphysics  Simulations:  Challenges  and  Opportunities,  D.  E.  Keyes,  L.  C.  Meinnes,  C. 
Woodward,  W.  Gropp,  E.  Myra,  M.  Pernice,  J.  Bell,  J.  Brown,  A.  Clo,  J.  Connors,  E.  Con- 
stantinescu,  D.  Estep,  K.  Evans,  C.  Farha^,  A.  Hakim,  G.  Hammond,  G.  Hansen,  J.  Hill, 

T.  Isaac,  X.  Jiao,  K.  Jordan,  D.  Kaushik,  E.  Kaxiras,  A.  Koniges,  K.  Lee,  A.  Lott,  Q.  Lu, 

J.  Magerlein,  R.  Maxwell,  M.  McCourt,  M.-Mehl,  R.  Pawlowski,  A.  Peters  Randles,  D. 
Reynolds,  B.  Riviere,  U.  Ruede,  T.  Scheibe,  J.  Shadid,  B.  Sheehan,  M.  Shephard,  A.  Siegel, 

B.  Smith,  X.  Tang,  C.  Wilson,  and  B.  Wohlmuth,  International  Journal  of  High  Performance 
Computing  Applications  (27),  2013. 
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•  Continuum  Modeling  and  Control  of  Large  Nonuniform  Wireless  Networks  via  Nonlinear 
Partial  Differential  Equations,  Y.  Zhang,  E.  Chong,  J.  Hannig,  and  D.  Estep,  Abstract  and 
Applied  Analysis  (16),  2013,  doi:10.1 155/2013/262581,  1-16 

•  A  posteriori  error  estimates  for  explicit  time  integration  methods,  J.  Collins,  D.  Estep  and  S. 
Tavener,  BIT  Numerical  Mathematics,  2012,  submitted 

•  Continuum  Limits  of  Markov  Chains  with  Application  to  Wireless  Network  Modeling,  Y. 
Zhang,  E.  Chong,  J.  Hannig,  and  D.  Estep,  IEEE  Access,  2013,  submitted 

•  A  posteriori  error  estimation  for  the  Lax-Wendroff finite  difference  scheme,  J.  B.  Collins,  D. 
Estep,  and  S.  Tavener,  Journal  of  Computational  and  Applied  Mathematics,  2013,  submitted 

•  Convergence  and  optimality  of  adaptive  methods  in  the  Finite  Element  Exterior  Calculus 
framework,  M.  Holst,  A.  Mihalik,  and  R.  Szypowski,  Submitted  for  publication. 

•  An  alternative  between  non-unique  and  negative  yamabe  solutions  to  the  conformal  formu¬ 
lation  of  the  einstein  constraint  equations,  M.  Holst  and  C.  Meier,  Submitted  for  publication. 

•  Non-uniqueness  of  solutions  to  the  conformal  formulation,  M.  Holst  and  C.  Meier,  Submitted 
for  publication. 

•  Efficient  computational  in  multiscale  geometric  modeling  for  biomolecular  complexes,  T.  Liao, 
Y.  Zhang,  P.  Kekenes-Huskey,  A.  Michailova,  M.  Holst,  and  J.  A.  McCammon,  Submitted 
for  publication. 

•  Multilevel  preconditioners  for  discontinuous  Galerkin  approximations  of  elliptic  problems 
with  jump  coefficients,  B.  Ayuso  de  Dios,  M.  Holst,  Y.  Zhu,  and  L.  Zikatanov,  Accepted  for 
publication  in  Math.  Comp. 

4.2.  Presentations  at  meetings,  conferences,  seminars 

The  following  presentations  were  made  during  March  27,  2009  -  September  1,  2009 

Burch:  Research  Seminar,  Sandia  National  Laboratory,  Albuquerque,  New  Mexico,  8/09 

Estep:  Computational  Science  and  Engineering  (CSE)  Annual  Research  Symposium,  University 
of  Illinois,  Urbana-Champaign,  Keynote  Speaker,  4/09 

Estep:  SIAM  Annual  Meeting,  Minisymposium  on  Predictive  Computational  of  Multiscale- 
Multiphysics  Applications,  invited  speaker,  7/09 

Estep:  Workshop  on  Simulating  the  Spatial-Temporal  Patterns  of  Anthropogenic  Climate  Change, 
Los  Alamos  Institute  for  Advanced  Studies,  Santa  Fe,  New  Mexico,  invited  speaker,  8/09 

Estep:  Colloquium,  Department  of  Mathematics,  University  of  Wyoming,  9/09 

Holst  25th  Pacific  Coast  Gravity  Meeting  (PCGM25),  Eugene,  Oregon,  4/09 

Holst:  5th  Annual  Structured  Integrators  Workshop,  Caltech,  Pasadena,  California,  Plenary  Speaker, 
5/09 

Holst:  FEniCS  2009  Workshop,  Oslo,  Norway,  Plenary  Speaker,  6/09 

Holst:  Numerische  Mathematik  50,  Munich,  Germany,  Plenary  Speaker,  6/09 

Holst:  Mathematical  and  Numerical  Geometric  Analysis  Workshop,  Frieburg,  Germany,  Plenary 
Speaker,  9/09 


Holst:  ICNAAM  Conference,  Crete,  Greece,  Minisymposium  Speaker,  9/09 
Serenesa:  CSME  Seminar  Series,  UC  San  Diego,  San  Diego,  California,  6/09 

The  following  presentations  were  made  during  September  2.  2009  -  September  1,  2010 
Burch:  ICMS  Workshop  on  Uncertainty  Quantification,  Edinburgh,  UK,  05/10 

Estep:  Workshop  on  Adaptive  and  Multilevel  Methods  for  Partial  Differential  Equations,  Univer¬ 
sity  of  California  San  Diego,  1 1/09 

Estep:  Seminar,  Lawrence  Livermore  National  Laboratory,  12/09 

Estep:  Colloquium,  Department  of  Atmospheric  Science,  Colorado  State  University,  1/10 

Estep:  Seminar,  University  of  Wisconsin,  2/10 

Estep:  Seminar,  Brown  University,  3/10 

Estep:  Seminar,  University  of  Chicago,  3/10 

Serenesa:  CCoM  Seminar  Series,  UC  San  Diego,  San  Diego,  California,  1 1/09 
Holst:^Pknary  Lecture,  Symposium  on  Mathematical  Systems  Biology,  UCI,  Irvine,  California, 

Holst:  Lecture,  26th  Pacific  Coast  Gravity  Meeting  (PCGM26),  San  Diego,  CA,  3/10 

Holst:  Plenary  Lecture,  Workshop  on  Unstructured  Mesh  Methods  in  Mathematical  Physics,  Jena, 
Germany,  8/10 

Holst:  Invited  Lecture,  Department  of  Mathematics,  Free  University  of  Berlin,  Berlin  Germany 
8/10 

Holst:  Invited  Lecture,  Department  of  Mathematics,  Technical  University  of  Berlin,  Berlin,  Ger¬ 
many,  8/10 

Holst:  Invited  Lecture,  Department  of  Mathematics,  Jacobs  University,  Brehmen,  Germany,  9/10 

The  following  presentations  were  made  during  September  2,  2010  -  September  1,2011 

Estep:  SIAM  Computational  Science  and  Engineering  Conference,  Minisymposia  on  Numerical 
Discretization  Error  Estimation  for  Uncertainty  Quantification,  Progress  in  Computational 
Methods  and  Software  for  Tightly-coupled  Multiphysics  Applications,  Numerical  Methods 
for  Stochastic  Computation  and  Uncertainty  Quantification,  Numerical  Challenges  in  Mi¬ 
crostructure  Modeling  for  Materials  Science,  Reno,  Nevada,  201 1 

Estep:  Seminar,  Lawrence  Livermore  National  Laboratory,  9/10 

Estep:  Seminar,  Purdue  University,  9/10  ''  ' 

Estep:  Seminar,  North  Carolina  State  University,  11/10 

Estep:  Seminar,  Lawrence  Livermore  National  Laboratory,  1/1 1 

Estep:  Seminar,  University  of  Southern  California,  3/1 1 

Estep:  Plenary  Talk,  ICiS  Workshop  on  Multiphysics  Simulations:  Challenges  and  Opportunities, 
Park  City,  Utah,  8/1 1 
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Holst:  Invited  Lecture,  Department  of  Mathematics,  Jacobs  University,  Bremen,  Germany,  9/10 

Holst:  Invited  Lecture,  Workshop  on  Latest  Trends  and  Developments  in  Computational  Tech¬ 
nology  and  Methods  for  Solids,  Structures,  Fluids  and  Fluid-Structure  Interaction,  La  Jolla 
CA,  9/10 

Holst:  Invited  ICES  Lecture,  University  of  Texas,  Austin,  TX,  2/1 1 
Holst:  Invited  CVS  Lecture,  University  of  Texas,  Austin,  TX,  2/1 1 

Holst:  Colloquium,  Department  of  Mathematics,  University  of  Wisconsin,  Madison,  WI,  4/1 1 

Holst:  Colloquium,  Department  of  Mathematics,  The  Penn  State  University,  State  College  PA 
4/11  /  6  . 

Holst:  Colloquium,  Department  of  Applied  Mathematics,  University  of  Washington,  Seattle,  WA, 


Holst:  Seminar,  Pacific  Northwest  National  Laboratory,  Richland,  WA,  5/1 1 

Holst:  Plenary  Lecture,  Workshop  on  Advances  and  Challenges  in  Computational  General  Rela¬ 
tivity,  Brown  University,  Providence,  RI,  5/1 1 

Holst.  Invited  Lecture,  Schnelle  Loser  fur  partielle  Differentialgleichungen,  Mathematisches 
Forschungsinstitut  Oberwolfach,  Oberwolfach,  Germany,  5/1 1 

The  following  presentations  were  made  during  September!,  2011  -  September  1,  2012 

Estep:  Invited  Lecture,  Uncertainty  Quantification  for  High-Performance  Computing  Workshop, 
Oak  Ridge  National  Laboratory,  5/12 

Estep:  Invited  Lecture,  6th  International  Conference  on  Automatic  Differentiation,  Fort  Collins 
CO,  7/12 

Estep:  Invited  Paper,  Joint  Statistical  Meetings,  8/12 

Estep:  Invited  Seminar,  University  of  Chicago,  9/1 1 

Estep:  Invited  Seminar,  Florida  State  University,  4/12 

Estep:  Invited  Seminar,  Colorado  School  of  Mines,  4/12 

Estep:  Invited  Colloquium,  Statistical  and  Applied  Mathematical  Sciences  Institute  (SAMSD 
4/12 

Holst:  Invited  Lecture,  Workshop  on  Geometric  Partial  Differential  Equations:  Theory,  Numer¬ 
ics  and  Appli-  cations,  Mathematisches  Forschungsinstitut  Oberwolfach,  Oberwolfach,  Ger¬ 
many,  11/11 

Holst:  Invited  Lecture,  JTO  Faculty  Fellowship  Lecture  (1  of  2),  Institute  for  Computational 
Engineering  and  Science  (ICES),  University  of  Texas,  Austin,  TX,  1 1/11 

Holst:  Invited  Lecture,  JTO  Faculty  Fellowship  Lecture  (2  of  2),  Institute  for  Computational 
Engineering  and  Science  (ICES),  University  of  Texas,  Austin,  TX,  1/12 

Holst:  Plenary  Lecture,  CSU  Research  Colloquium,  Physics  at  CSU:  Neutrinos  to  Nano  Science, 
Colorado  State  University,  Fort  Collins,  CO,  3/12 

Holst:  Plenary  Lecture,  21st  International  Conference  on  Domain  Decomposition  Methods,  Rennes 
Frances,  6/12  ,  ,  . 
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The  following  presentations  were  made  during  September  2,  2012  -  April  5,  2013 

Pollock:  Center  for  Computational  Mathematics  Seminar,  UCSD,  San  Diego,  CA,  1/13. 

Pollock:  Joint  MAA-AMS  Mathematics  Meetings,  San  Diego,  CA,  1/13. 

Pollock:  Numerical  analysis  seminar,  Texas  A&M  University,  College  Station,  TX,  4/13. 

Pollock:  CSME  Seminar,  UCSD,  San  Diego,  CA,  4/13. 

Pollock:  Minisymposium  Lecture,  SIAM  Annual  Meeting,  San  Diego,  CA,  7/13. 

4.3.  Websites 

Research  results  and  software  are  presented  at 

•  http://www.stat.colostate.edu/~estep/ 

•  http://ccom.ucsd.edu/~mholst/ 

4.4.  Technologies  and  techniques 

Over  the  last  several  years,  our  DTRA-supported  research  team  has  led  the  development  of 
the  Finite  Element  ToolKit,  which  is  an  opensource  finite  element  modeling  toolkit  designed  for 
the  simulation  of  coupled  multiphysics  problems  with  multiscale  phenomena.  The  software  has 
been  designed  and  developed  collaboratively  by  both  Holst  and  Estep,  and  consists  of  a  collection 
of  object-oriented  class  libraries  written  in  C,  C++,  Objective  C,  and  Python.  There  is  also  a 
MATLAB/Octave-based  prototyping  tool  (FETKLab),  the  development  of  which  has  been  done 
by  both  Estep  and  Holst,  as  well  as  several  of  their  graduate  students.  FETK  (and  FETKLab)  are 
designed  to  adaptive  discretize  and  solve  coupled  reaction-diffusion  systems,  and  is  based  around 
state-of-the-art  algorithms  for  simplex  mesh  generation,  error  estimation,  mesh  refinement,  finite 
element  discretization,  iterative  nonlinear  and  optimization  techniques,  and  fast  multilevel  and 
domain  decomposition-based  linear  solvers  and  preconditions.  Many  of  the  algorithms  developed 
in  our  research  articles  as  described  in  this  report  have  been  prototyped,  implemented,  and  applied 
to  applications  in  conjunction  with  physical  scientists  using  FETK.  The  entire  FETK  source  tree 
was  released  in  June  2010  on  the  FETK.org  website,  as  a  major  milestone  of  this  DTRA  project. 
A  substantial  extention  to  both  FETK  and  FETKLab  was  completed  in  Spring  2013  that  added 
general  Lagrange-type  elements  for  either  primal  or  dual  problems,  and  this  new  capability  has 
been  exploited  in  a  number  of  our  recent  articles. 

In  addition,  we  continue  development  on  GAASP  (Globally  Accurate,  Adaptive  Sensitivity 
analysis  Package)  to  extend  its  capabilities  for  both  forward  and  inverse  stochastic  sensitivity  anal¬ 
ysis  of  differential  equations. 

5.  Impact 

5.1.  Impact  on  the  principal  disciplines  of  the  project 

The  numerical  solution  of  multiscale,  multiphysics  models  on  complex  domains  along  with 
the  development  of  tools  for  predictive  science  and  uncertainty  quantification  is  one  of  the  grand 
challenges  facing  the  mathematical  sciences  at  present.  Such  problems  present  a  very  complex  pic¬ 
ture  in  terms  of  stability  and  important  behaviors  interacting  across  a  wide  range  of  scales,  which 
makes  the  straightforward  use  of  classical  numerical  methods  and  analyses  extremely  problematic, 
if  not  impossible.  Classic  approaches  were  developed  in  the  context  of  models  involving  single 
physics  phenomena  operating  at  a  narrow  range  of  scales.  While  building  on  classic  approaches, 
the  research  in  this  project  contributes  at  a  fundamental  theoretical  level  by  laying  the  foundation 
for  reliably  accurate  and  efficient  numerical  solution  based  on  a  posteriori  error  analysis  that  ac¬ 
counts  for  the  numerical  complexities  involved  with  simulating  such  systems.  This  is  achieved 
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by  combining  extremely  sophisticated  mathematics  in  analysis  and  geometry  with  cutting  edge 
numerical  methodology. 

The  impact  of  the  research  related  to  this  project  is  widespread,  as  can  be  seen  in  the  greatly 
increasing  levels  of  activity  around  the  world  on  such  problems.  This  is  also  evidenced  by  the 
number  of  invitations  to  speak,  the  number  of  funded  interdisciplinary  projects  including  a  recent 
award  of  an  extremely  prestigious  National  Science  Foundation  Focused  Research  Group  (FRG) 
award  to  Estep  and  Holst,  the  citation  record  (Estep’s  h-index  is  15  and  Holst’s  h-index  is  20), 
and  the  high  level  of  the  involvement  of  the  Pi’s  in  research  environment  through  panels,  reports, 
editing,  and  so  on. 

5.2.  Impact  on  other  disciplines 

Developing  reliable  and  accurate  tools  for  carrying  out  predictive  science  and  engineering  for 
multiscale,  multiphysics  systems  on  complex  domains  and  conducting  uncertainty  quantification 
in  simulated  results  is  the  major  problem  of  computational  science  and  engineering  at  present. 
Addressing  this  challenge  requires  fundamental  research  in  the  mathematical  sciences.  This  project 
is  aimed  at  addressing  a  number  of  key  research  problems  involved  with  simulating  multiphysics 
systems.  Along  with  theory,  the  Pis  systematically  implement  the  results  into  public  software, 
and,  along  with  their  collaborators,  use  the  software  to  tackle  .scientific  and  engineering  research 
problems.  This  yields  a  direct  transfer  of  the  theoretical  mathematical  developments  and  software 
implementations  to  the  application  domain. 

This  is  evidenced  by  the  large  number  of  interdisciplinary  collaborations  of  the  Pis  and  the 
substantial  volume  of  interactions  with  Department  of  Energy  laboratories  and  industry.  Details 
are  provided  below. 

5.3.  Impact  in  the  profession 

5.4.  Honors  and  awards 

Estep  was  appointed  (founding)  Co-Editor  in  Chief  of  the  SIAM  /  ASA  Journal  on  Uncertainty 
Quantification 

Estep  won  the  University  Scholarship  Impact  Award,  Colorado  State  University,  201 1 

Estep  was  appointed  University  Interdisciplinary  Research  Scholar,  Colorado  State  Universitv  in 
2009  ^ 

Estep  received  the  Oliver  P.  Pennock  Distinguished  Service  Award,  Colorado  State  Universitv  in 
2009  ^ 

Estep  was  appointed  Editor  in  Chief,  SIAM  Book  Series  on  Computational  Science  and  Engi¬ 
neering,  2009  -  2014 

Holst  received  the  CSU  Distinguished  Alumnus  Award,  2009 

Holst  was  appointed  the  Chancellor’s  Associates  Endowed  Chair  in  Mathematics  and  Physics  at 
UC  San  Diego  in  2012 

5.5.  Impact  on  the  professional  research  community 

Estep  served  as  one  of  the  Moderators  for  the  SAMSI  National  SIAM  and  ASA  Town  Hall  Meet¬ 
ing  on  Uncertainty  Quantification,  2010 

Estep  served  as  the  Co-Organi/er  and  first  Chair,  SIAM  Activity  Group  on  Uncertainty  Quantifi¬ 
cation,  2010 

Estep  served  as  a  Program  Leader  for  the  SAMSI  Program  on  Uncertainty  Quantification,  2011- 
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Estep  served  as  co-chair  of  the  first  SIAM/ASA/USACM  Conference  on  Uncertainty  Quantifica¬ 
tion  (April,  2012) 

Estep  along  with  J.  Berger  (Duke)  and  M.  Gunzburger  (FSU)  proposed  a  new  Journal  on  Uncer¬ 
tainty  Quantification  to  be  jointly  published  by  the  ASA  and  SIAM 

Estep  serves  on  the  Advisory  Board  for  the  Center  for  Advanced  Modeling  and  Simulation,  Idaho 
National  Laboratory,  2009  -  2012 

Estep  serves  on  the  Governing  Board  of  the  Statistical  and  Applied  Mathematical  Sciences  Insti¬ 
tute  (SAMSI),  2009-2016 

Estep  served  on  the  National  Science  Foundation  Office  of  Cyberinfrastructure  Grand  Challenges 
Communities  Task  Force,  2009-2010  (co-author  of  final  recommendation  report) 

Estep  served  as  Breakout  Lead  and  Report  co-author.  Uncertainty  Quantification  and  Stochastic 
Systems,  Department  of  Energy  Cross-Cutting  Technologies  for  Computing  at  the  Exascale 
2010 

Estep  was  an  invited  participant  in  the  Fusion  Simulation  Program  Definition  Workshop,  201 1 

Estep  serves  on  the  American  Mathematical  Society  Simmons  Travel  Grants  Committee  2011- 
2014 

Estep  serves  as  Moderator,  Mathematics  in  the  Geosciences  Workshop,  Northwestern  University, 

Estep  was  co-author  of  Multiphysics  Simulations:  Challenges  and  Opportunities,  Tech.  Report 
ANL/MCS-TM-321,  Argonne  National  Laboratory,  201 1 

Estep  was  co-author  of  Fostering  Interactions  Between  the  Geosciences  and  Mathematics,  Statis¬ 
tics,  and  Computer  Science,  Technical  Report  TR-20 12-02,  Department  of  Computer  Sci¬ 
ence,  University  of  Chicago,  2012 

Holst  serves  on  the  Executive  Committee  for  the  San  Diego  Supercomputer  Center  (SDSC),  2007- 
present 

Holst  is  a  Co-Organizer  (with  R.  Bank)  of  20th  International  Conference  on  Domain  Decompo¬ 
sition  (DD20),  February  2011. 

Holst  is  the  Primary  Organizer  (with  J.  Hameed):  Numerical  Methods  for  Implicit  Models  in 
Biomolecular  Systems,  SIAM  CS&E  Conference  Minisymposium,  March  201 1 

Holst  is  the  Primary  Organizer  (with  A.  Demlow,'  A.  Gillette,  Y.  Zhu):  Workshop  on  Exploit¬ 
ing  Geometry  in  the  Development  of  Numerical  Methods  for  Partial  Differential  Equations, 
UCSD  Workshop,  San  Diego,  November  2011. 

Holst  is  the  Primary  Organizer  (with  A.  Demlow,  R.  Szypowski):  Exploiting  Geometry  in  the 
Development  of  Numerical  Methods  for  Partial  Differential  Equations,  SIAM  Analysis  of 
PDE  Conference  Minisymposium,  November  2011. 

Holst  is  the  Primary  Organizer  (with  D.  Arnold,  A.  Gillette):  AMS  Joint  Meeting  FEEC  Min¬ 
isymposium,  on  New  Developments  in  the  Finite  Element  Exterior  Calculus,  January  2013. 

Holst  is  the  Primary  Organizer  (with  A.  Gillette,  R.  Szypowski):  Workshop  on  Exploiting  Geom¬ 
etry  in  the  Development  of  Numerical  Methods  for  Partial  Differential  Equations  II,  UCSD 
Workshop,  San  Diego,  January  2013. 

Holst  and  Estep  regularly  serve  on  Grant  Review  Panels  for  NSF  and  DOE,  2004-present 
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5.6.  Professional  editorial  appointments 

Estep:  CO  Editor  in  Chief  (founding),  SIAM  /  ASA  Journal  on  Uncertainty  Quantification 

Estep:  Editor  in  Chief,  SIAM  Book  Series  on  Computational  Science  and  Engineering,  2009  - 
2014 

I  I 

Estep:  Associate  Editor,  SIAM  Journal  on  Numerical  Analysis,  2005-201 1 

Estep:  Associate  Editor,  International  Journal  for  Uncertainty  Quantification,  2010- 

Estep:  Associate  Editor,  Multiphysics  Modeling  Book  Series,  A.  A.  Balkema  Publishing,  CRC 
Press,  2010- 

Estep:  Associate  Editor,  Journal  of  Applied  Mathematics  and  Computing,  2008-2013 
Holst:  Associate  Editor,  Numerische  Mathematik,  2008-present 

Holst:  Associate  Editor,  SIAM  Book  Series  on  Computational  Science  and  Engineering  2009- 
2014 

5.7.  Impact  on  technology  transfer 

The  Pis  maintain  a  very  substantial  interdisciplinary  collaboration  activity  with  scientists  and 
engineers  in  universities.  Department  of  Energy  laboratories,  and  industry.  These  collaborations 
lead  to  direct  injection  of  research  ideas  into  practical  use. 

5.8.  Consulting  and  collaborative  activities 

In  this  section,  we  report  currently  funded  projects  that  involve  substantial  interdisciplinary 
collaborations  and  transfer  of  research  results  related  to  this  project  into  applications. 

Estep  IS  CO-PI  on  the  project  Framework  Application  for  Core-Edge  Transport  Simulations  (FACETS) 
funded  by  the  Office  of  Advanced  Scientific  Computing  Research  and  Office  of  Fusion  En¬ 
ergy  Sciences,  Department  of  Energy,  2007-12.  Collaborators  include:  R.  H.  Cohen,  L.  Di- 
achin,  and  T.  Eppcrly  at  Lawrence  Livermore  National  Laboratory;  J.  Larson  and  L.  Mclnnes 
at  Argonne  National  Laboratory;  M.  R.  Fahey  and  J.  Cobb  at  Oak  Ridge  National  Laboratory. 
Subject  is  development  and  analysis  of  numerical  solution  methods  for  coupled  core-edge 
fusion  simulations. 

Estep  is  PI  on  the  project  Collaborative  Proposal:  Transforming  How  Climate  System  Models 
are  Used:  A  Global,  Multi- Resolution  Approach  to  Regional  Ocean  Modeling  funded  by  the 
Department  of  Energy,  2009-1 1 .  Collaborators  include  Todd  Ringler  at  Los  Alamos  National 
Laboratory.  Subject  is  development  and  analysis  of  numerical  methods  for  multiscale  ocean 
models. 

Estep  is  PI  on  the  project  Adjoint-based  methods  for  uncertainty  quantification  funded  by  the 
Lawrence  Livermore  National  Laboratory,  2010-13.  Collaborators  are  Carol  Woodward  and 
Jeff  Hittinger  at  Lawrence  Livermore  National  Laboratory.  Duties  include  (1)  pursue  develop 
a  posteriori  error  estimates  for  hyperbolic  problems  including  shock  behavior  and  (2)  consult 
on  uncertainty  and  error  quantification  with 'laboratory  personnel 

I  ] 

Estep  is  CO-PI  on  the  project  The  Inverse  Problem  for  Estimation  of  Structure  of  Biological  Macro¬ 
molecules  from  Small-Angle  X-Ray  Scattering  funded  by  the  National  Institutes  of  Health, 
2010-2014.  Collaborators  include  Jay  Breidt  (Statistics,  CSU)  and  Karolin  Luger  (Biochem¬ 
istry,  CSU).  Subject  is  determining  the  structure  of  biological  macromolecules  using  small 
angle  x-ray  scattering  data. 
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Estep  is  PI  on  the  project  Enabling  Predictive  Simulation  and  UQ  of  Complex  Multiphysics 
PDE  Systems  by  the  Development  of  Goal-Oriented  Variational  Sensitivity  Analysis  and 
a-Posteriori  Error  Estimation  Methods  funded  by  the  Department  of  Energy,  2010-2013. 
Collaborators  include  John  Shadid  (Sandia  Nat.  Lab.)  and  Victor  Ginting  (U.  Wyom.)' 
Subject  is  developing  a  posteriori  error  estimates  for  solutions  of  reacting  flow  and  fusion 
reaction  models. 

Estep  is  CO-PI  on  the  project  Collaborative  Research:  A  posteriori  error  analysis  and  adaptiv¬ 
ity  for  discontinuous  interface  problems  funded  by  the  National  Science  Foundation,  2010- 
20n.  Collaborator  is  Simon  Tavener  (CSU).  Purpose  is  developing  and  analyzing  conser¬ 
vative  solution  methods  for  elliptic  problems  with  coefficients  that  are  discontinuous  across 
complex  interfaces. 

Estep  is  PI  on  the  CSU  Subcontract  from  Multiscale  Design  Systems,  LLC  supported  by  an 
Air  Force  SBIR  Phase  II  grant.  Collaborators  are  Simon  Tavener  (CSU)  and  Jacob  Fish 
(Columbia  Uni.)  in  2011.  Purpose  is  developing  fast  methods  for  UQ  for  multiscale  models 
of  polymers  in  stressed  environments. 

Estep  is  PI  on  the  project  Uncertainty  Analysis  for  Multiscale  Models  of  Nuclear  Fuel  Perfor¬ 
mance  supported  by  the  Idaho  National  Laboratory  from  201 1-2014.  Collaborators  are  Si¬ 
mon  Tavener  (CSU)  and  Michael  Pernice  (Idaho  Nat.  Lab.).  Purpose  is  UQ  for  multiscale 
models  of  nuclear  fuel. 

Estep  is  PI  on  the  project  11-2031:  Multiscale  modeling  and  uncertainty  quantification  for 
nuclear  fuel  performance.  Nuclear  Energy  University  Programs,  Department  of  Energy, 
2011-14.  Collaborators  are  Simon  Tavener  (CSU),  Michael  Pernice  (INL),  Peter  Polyakov 
(Wyoming),  Dongbin  Xiu  (Purdue),  Anter  el  Azab  (Purdue) 

Estep  is  a  co-PI  on  the  project  Data-Driven  Inverse  Sensitivity  Analysis  for  Predictive  Coastal 
Ocean  Modeling,  Computational  and  Data-Enabled  Science  and  Engineering  in  Mathemati¬ 
cal  and  Statistical  Sciences  (CDS&E-MSS),  National  Science  Foundation,  2012-15.  Collab¬ 
orators  are  Troy  Butler  (CSU),  Clint  Dawson  (U.  Texas  at  Austin),  and  Joannes  Westerink 
(Notre  Dame) 

Estep  and  Holst  are  co-PIs  on  the  project  ERG:  Error  Quantification  and  Control  for  Gravita¬ 
tional  Waveform  Simulation  funded  by  the  National  Science  Foundation,  2011-2014.  The 
Project  is  concerned  with  estimating  the  error  in  computed  wave  forms  obtained  from  LIGO 
data.  '  • 

Holst  is  Pis  on  the  project  hRG:  Analysis  of  the  Einstein  Constraint  Equations  funded  by  the 
National  Science  Foundation,  2013-2016.  The  Project  is  concerned  with  further  extending 
the  solution  theory  for  the  Einstein  constraint  equations. 

Holst  is  PI  on  the  project  MRI:  Acquisition  of  a  Parallel  Computing  and  Visualization  Facility  to 
Enable  Integrated  Research  and  Training  in  Modern  Computational  Science,  Mathematics 
and  Engineering  funded  by  National  Science  Foundation,  2008-2011.  Collaborators  include 
Randolph  Bank  (UCSD  Mathematics),  Scott  Baden  (UCSD  Computer  Science),  and  John 
Weare  (UCSD  Chemistry).  The  subject  is  the  design  and  construction  of  a  state-of-the-art 
parallel  computing  system  with  an  excess  of  1000  compute  nodes,  Infiniband  high-speed 
network  fabric,  parallel  filesystems,  LCD  vizualization  walls,  housed  in  a  modern  server 
room  with  raised  floor  and  forced  chilled  air. 

Holst  is  PI  on  the  project  Adaptive  Methods  and  Finite  Element  Exterior  Calculus  for  Nonlinear 
Geometric  PDE,  funded  by  National  Science  Foundation,  2012-2015.  Co-PI  is  former  stu¬ 
dent  and  postdoc  Ryan  Szypowski,  now  an  assistant  professor  in  mathematics  at  Cal  Poly 
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Pomona.  The  subject  is  the  design  and  analysis  of  adaptive  methods  for  use  with  the  finite 
element  exterior  calculus. 

Holst  is  Co-PI  on  the  project  Adaptive  Radiotherapy  Based  on  High  Performance  Computing 
funded  by  the  Department  of  Energy,  Lawrence  Livermore  National  Laboratory,  and  the  Uni¬ 
versity  of  California,  2009-2012.  Collaborators  include  Steve  Jian  (UCSD  Medical  School) 
A.  Majumdar  (SDSC),  and  D.J.  Choi  (SDSC).  The  subject  is  realtime  solution  of  coupled 
reaction-diffusion  systems  and  the  Boltzmann  transport  equation  using  a  combination  of  par¬ 
allel  algorithms  for  partial  differential  equations,  high-speed  communication  networks,  and 
cluster  computers. 

Holst  is  Co-PI  on  the  project  Scalable  Adaptive  Multilevel  Solvers  for  Multiphysics  Problems, 
funded  by  the  Department  of  Energy.  The  subject  is  the  design  and  analysis  of  determinstic 
algorithms  for  use  in  physical  simulation  based  on  multilevel  technologies. 

Holst  is  Co-PI  on  the  project  Applications  of  Quantum  Computing  in  Aerospace  Science  and 
Engineering,  funded  by  the  AirForce  Office  of  Scientific  Research.  The  subject  is  the  design 
and  analysis  of  quantum  algorithms  for  use  in  physical  simulation. 

Holst  is  Co-PI  and  Core  lA  lead  on  the  project  National  Biomedical  Computation  Resource 
(NBCR)  funded  by  the  National  Institutes  of  Health,  2009-2014.  Collaborators  include  An¬ 
drew  McCammon  (UCSD  Chemistry),  Andrew  McCulloch  (UCSD  Bioengineering),  Mark 
Elhsman  (UCSD  Medical  School),  and  Peter  Arzberger  (SDSC).  The  subject  is  multiscale 
modeling  frameworks  and  adaptive  finite  clement  methods  for  complex  multiscale  and  mul¬ 
tiphysics  problems  arising  in  biomedical  science. 

Holst  is  Senior  Scientist  and  founding  member  of  the  NSF  Physics  Frontier  Center  for  Theoret¬ 
ical  Biological  Physics  (CTBP),  funded  by  the  National  Science  Foundation.  Collaborators 
include  Jose’  Onuchic  (UCSD  Physics),  Andrew  McCammon  (UCSD  Chemistry),  and  Andy 
Kummel  (UCSD  Chemistry).  The  subject  is  multiscale  modeling  frameworks  and  adaptive 
finite  element  methods  for  complex  multiscale  and  multiphysics  problems  arising  in  bio¬ 
physics. 

Transitions  to  technology  applications 
We  report  on  current  interactions  with  industry. 

Estep  was  a  Co-Principal  Investigator  in  the  Tech  X,  Inc.  project  Framework  Application  for 
Core-Edge  Transport  Simulations  (FACETS),  funded  by  the  Office  of  Advanced  Scientific 
Computing  Research  and  Office  of  Fusion  Energy  Sciences,  Department  of  Energy.  Estep’s 
responsibilities  include  development  and  analysis  of  numerical  solution  methods  for  coupled 
core-edge  fusion  simulations.  Algorithms  developed  in  this  program  will  be  implemented 
into  the  FACETS  high  performance  framework. 

Estep  was  a  subcontract  in  Phase  II  project  for  Multiscate  Design  Systems,  LLC  (Principal  Officer; 
Jacob  Fish,  Rensselaer  Polytechnic  Institute)  for  the  Air  Force  SBIR/STTR  program.  Es¬ 
tep  s  responsibilities  include  development  of  niultiscale  operator  decomposition  numerical 
methods  and  numerical  methods  for  error  estimation,  uncertainty  quantification  and  inverse 
problems  for  parameter  identification  for  multiscale  multiphysics  models  of  hygro-thermo- 
mechano-oxidation-fatigue  in  polymer  matrix  composites  used  in  aircraft  applications.  Al- 
gonthms  developed  in  this  program  will  be  implemented  into  the  Multiscale  Design  Sys¬ 
tem  for  Continuum  (MDS-C)  and  the  Multiscale  Design  System  for  Discrete  or  atomistic 
medium  (MDS-D)  software  packages. 
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Holst  is  collaborating  with  Eric  Bylaska  at  Pacific  Northwest  National  Laboratory  on  the  incorpo¬ 
ration  of  the  Finite  Element  Toolkit  (FETK,  developed  and  maintained  by  the  Holst  Group) 
into  several  density  functional  modeling  packages  based  at  PNNL. 


24 


BLOCKWISE  ADAPTIVITY  FOR  TIME  DEPENDENT  PROBLEMS 
BASED  ON  COARSE  SCALE  ADJOINT  SOLUTIONS 

V,  CAREY  D.  ESTEP  t,  A.  JOHANSSON  *,  M.  LARSON  5,  and  S.  TAVENER  ’ 

Abstract.  We  describe  and  test  an  adaptive  algorithm  for  evolution  problems  that  employs  a 
sequence  of  “blocks”  consisting  of  fixed,  though  non-uniform,  space  meshes.  This  approach  offers 
the  advantages  of  adaptive  mesh  refinement  but  with  redueeri  overhead  costs  associated  with  load 
balancing,  re-meshing,  matrix  reassembly,  and  the  solution  of  adjoint  problems  used  to  estimate 
discretization  error  and  the  cfTccts  of  mesh  changes.  A  major  issue  with  a  block-adaptive  approach 
is  determining  block  discretizations  from  coarse  scale  solution  information  that  achieve  the  desired 
accuracy.  We  de.scribe  several  strategies  to  achieve  this  goal  using  adjoint-based  a  posteriori  error 
estimates  and  we  demonstrate  the  behavior  of  the  proposed  algorithms  as  well  as  several  technical 
issues  in  a  set  of  examples. 

Key  words,  a  posteriori  error  analysis,  adaptive  error  control,  adaptive  mesh  refinement, 
adjoint  problem,  discontinuous  Galerkin  method,  duality,  generalized  Green’s  function,  goal  oriented 
error  estimates,  residual,  variational  analysis 

AMS  subject  classifications.  65N15,  65N30,  65N50 

1.  Introduction.  We  describe  and  test  an  adaptive  algorithm  for  evolution 
problems  that  we  call  “blockwise  adaptivity”.  This  approach  employs  a  sequence 
of  “blocks”  consisting  of  fixed,  though  non-uniform,  space  meshes,  and  is  motivated 
by  considerations  of  efficiency  and  accuracy.  We  balance  the  goal  of  achieving  de¬ 
sired  accuracy  using  discretizations  with  relatively  few  degrees  of  freedom  against  the 
computational  costs  iissociated  with  load  balancing,  re-meshing,  matrix  reassembly 
and  in  particular  the  cost  of  error  estimation.  A  block  adaptive  strategy  reduces  the 
number  of  mesh  changes  that  must  be  treated,  which  reduces  the  amount  of  com¬ 
putational  time  spent  on  re-meshing,  assembly,  and  loaxi  balancing,  and  makes  the 
problem  of  quantifying  the  effects  of  mesh  changes  on  accuracy  computationally  fea¬ 
sible.  A  block  adaptive  strategy  also  provides  a  natural  coarse  scale  discretization 
on  which  to  solve  the  adjoint  problem  used  to  compute  global  a  posteriori  error  esti¬ 
mates.  This  reduces  the  twin  computational  difficulties  of  storing  a  fine  scale  forward 
solution  in  order  to  form  the  adjoint  problem  and  solving  the  adjoint  problem  on  that 
fine  scale  discretization.  However,  a  major  issue  with  a  block-adaptive  approach  is 
determining  block  discretizations  from  coarse  scale  solution  information  that  achieve 
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the  desired  aecuracy  and  efiicioncy.  We  describe  several  strategics  to  achieve  this  goal 
using  adjoint-based  a  posteriori  error  estimates. 

To  focus  the  discussion,  we  consider  a  reaction-diffusion  equation  for  the  solution 
u  on  an  interval  [0,TJ, 

u  -  V  •  (£(a;,f)Vu)  =  f(u,x,t),  {x,t)  6  ff  x  (0,T], 

<  u{x,  t)  =  0,  (x,  t)  G  on  X  (0,  T),  (1.1) 

u(x,0)  =  Uo(x),  X  e  fi, 

where  fl  is  a  convex  polygonal  domain  in  with  boundai'y  dCl,  ii  denotes  the  partial 
derivative  of  u  with  respect  to  time,  and  there  is  a  constant  c  >  0  such  that 

f(x,  <)  >  f,  X  c  il,  t  >  0. 

We  also  assume  that  c  and  f  have  smooth  second  derivatives.  The  algorithms  in  this 
paper  generalize  to  problems  with  different  boundary  conditions,  convection,  nonlinear 
diffusion  coefficients,  as  well  as  systems,  sec  [17,  15], 

In  terms  of  adaptive  mesh  refinement,  the  interesting  situation  is  a  solution  of 
(1.1)  that  exhibits  "regionalized”  behavior  in  space  and  time.  Considerations  of  effi¬ 
ciency  suggests  that  time  steps  and  space  meshes  should  be  locally  refined  to  match 
the  regional  behavior,  see  the  plot  on  the  left  in  Fig.  1.1.  Classic  adaptive  mesh  re¬ 
finement  can  be  described  as  a  constrained  optimization  problem,  e.g.,  determine  a 
discretization  using  the  fewest  degrees  of  freedom  that  yields  a  solution  satisfying  a 
given  error  criterion.  In  general,  it  is  impossible  to  determine  a  clo.sed-form  solution 
of  this  optimization  problem.  An  adaptive  algorithm  is  an  iterative  procedure  for 
determining  a  nearly  optimal  solution. 


Eio.  1.1.  The  evolution  of  a  traveling  front  solution.  Left:  A  computation  using  space  meshes 
chosen  by  a  standard  adaptive  strategy  to  control  the  spatial  residual  error  at  each  time  step.  This 
entails  re-meshing,  re-asscmbly,  load  balancing,  and  projecting  the  solution  on  a  new  mesh  at  each 
step.  Right:  The  uniform  mesh  that  is  required  to  achieve  the  same  control  over  the  residual.  The 
computation  is  assembled  and  load  balanced  only  once. 

We  present  a  generic  adaptive  algorithm  in  Algorithm  1.1.  An  adaptive  compu¬ 
tation  is  generally  started  with  an  initial  coarse  mesh.  The  adaptive  algorithm  is  then 
applied  “real-time”  as  the  integration  proceeds  so  as  to  generate  a  new  space  mesh 
for  each  new  time  step,  where  the  new  space  mesh  is  based  on  (or  adapted  from)  the 
mesh  for  the  current  time  step.  In  practice,  the  remeshing  may  be  applied  on  intervals 
of  a  small  number  of  steps. 

While  adaptive  mesh  refinement  is  appealing  on  an  intuitional  level,  there  are 
serious  issues  facing  its  use  for  evolution  problems  including  the  following. 
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Algorithm  1.1  Generic  Adaptive  Algorithm  for  an  Evolution  Problem 
1:  Choose  an  initial  coarse  mesh  and  time  step 
2:  while  the  final  time  has  not  been  reached  do 

3:  Compute  a  numerical  solution  using  the  current  time  step  and  space  mesh 

4:  Estimate  the  error  of  the  computed  solution 

5:  while  the  error  estimate  is  too  large  do 

6:  Estimate  local  error  contributions  and  adapt  in  space 

7:  Estimate  local  error  contributions  and  adapt  in  time 

8;  Compute  a  numerical  solution  using  the  new  time  step  and  space  mesh 

9:  Estimate  the  error  of  the  computed  solution 

10:  end  while 

11:  Increment  time  by  the  accepted  time  step 

12:  end  while 


1.  Accuracy  Each  spatial  mesh  change  requires  a  projection  of  the  numerical 
sohition  onto  the  new  mesh,  and  this  can  affect  accuracy.  In  fact,  this  can 
destroy  convergence  altogether,  see  [8J. 

2.  Overhead  Costs  Changing  the  spatial  discretization  requires  generating  a 
new  mesh  and  reassembly  of  matrices.  Significant  mesh  changes  require  a 
redistribution  of  unknowns  among  the  processors  to  achieve  load  balancing. 
All  of  these  tasks  are  computational!}'  intensive. 

3.  Coarsening  Un-refinemcnt  or  coarsening  of  a  mesh  involves  loss  of  informa¬ 
tion  about  a  numerical  solution  that  cannot  be  recovered.  Currently,  there  is 
no  theory  for  coarsening  that  guarantees  that  there  is  no  loss  of  accuracy. 

4.  Global  Error  Estimation  Efficient  adaptive  mesh  refinement  requires  fic- 
curatc  error  estimates  of  the  true,  global  error,  but  cancelation  of  errors  over 
both  space  and  time  makes  choosing  adapted  meshes  problematic. 

Using  a  fixed  spatial  mesh  eliminates  the  first  three  issues.  But,  the  scale  required  of 
the  mesh  is  determined  by  the  finest  scale  required  in  any  region  where  discretization 
impacts  global  accuracy,  sec  Pig.  1.1.  This  necessarily  increases  computational  time 
and  solver  costs  and  memory  limits  may  make  it  impossible  to  use  the  necessary 
uniform  mesh. 

In  this  paper,  we  propose  a  “blockwise”  adaptive  algorithm  that  employs  nonuni¬ 
form  meshes  that  remain  fixed  for  discrete  period  of  times,  or  “blocks”,  sc'c  Fig.  1.2. 
With  the  proper  implementation,  this  strategy  addresses  the  following  key  issues. 

1.  Accuracy  The  projections  onto  new  meshes  occur  at  a  relatively  small  set 
of  discrete  times.  We  use  a  posteriori  error  estimates  to  predict  the  effect  of 
the  projections  and  choose  overlaps  in  the  meshes  to  reduce  the  error  induced 
by  the  mesh  changes. 

2.  Overhead  Costs  Re-meshing,  assembly,  and  load  balancing  are  required 
only  at  the  discrete  times  demarcating  blocks. 

3.  Coarsening  There  is  no  coarsening  of  a  given  mesh  in  the  indicated  strategy. 
Mesh  changes  are  handled  purely  as  projections  between  different  meshes. 

The  idea  of  rc-meshing  only  after  a  fixed  number  of  steps  is  by  no  means  new. 
However,  this  strategy  depends  critically  upon  choosing  suitable  block  discretizations, 
and  thus,  ultimately,  on  accurately  predicting  the  behavior  of  the  solution.  The  choice 
of  block  discretizations  is  a  difficult  issue  that  requires  balancing  the  inefficiency  of 
using  a  fixed  spatial  mesh  inside  each  block  against  the  gain  in  acc-uracy  achieved 
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Mesh 


Fig.  1.2.  The  evolution  of  a  solution  with  a  traveling  front  computed  zising  blockwise  adaptivity 
with  two  blocks.  On  each  block,  the  space  mesh  is  chosen  to  maintain  the  same  level  of  control 
over  the  local  residual  as  is  achieved  in  the  computation  shoum  in  Fig.  l.l.  In  addition,  there  is  a 
sufficient  degree  of  overlap  between  the  two  meshes  (the  lightly-shaded  mesh  region)  to  insure  there 
is  no  loss  of  accuracy  in  projecting  the  solution  between  the  two  meshes.  Re-meshing,  assembly,  and 
load  balancing  is  only  required  twice,  once  for  each  block. 


by  limiting  projections  between  different  meshes  and  the  decrease  in  computational 
cost  due  to  limiting  the  number  of  times  at  which  re-meshing,  re-assembly,  and  load 
balancing  is  required.  This  is  partly  a  computer  science  problem  of  distributing  avail¬ 
able  resources,  e.g.,  memory  and  compute  cycles,  efficiently,  and  partly  a  numerical 
analysis  problem,  e.g.,  determining  meshes  for  each  block  and  projections  between 
blocks. 

In  this  paper,  we  focus  on  the  problem  of  determining  blocks,  e.g.,  the  length 
of  times  for  each  block,  the  meshes  for  each  block  that  maintain  accuracy  in  the 
desired  information,  and  suitable  overlap  meshes  for  transitions  between  blocks  from 
the  coarse-scale  adjoint  solutions.  The  solutions  of  these  problems  require  accurate 
estimates  of  the  error  in  a  specific  quantity  of  interest.  We  use  a  computable  a 
posteriori  error  estimate  that  yields  robustly  accurate  estimates  of  the  error  in  a 
specified  quantity  of  interest  in  terms  of  a  sum  of  space-time  clement  contributions, 
see  [9,  10,  17,  15,  3,  20].  The  a  posteriori  error  estimates  are  based  on  duality, 
adjoint  problems,  and  variational  analysis.  Accurate  error  estimates  are  obtained 
by  numerically  solving  the  linear  adjoint  problem  related  to  the  desired  quantity  of 
interest. 

Solving  adjoint  problems  offers  computational  challenges  such  as  the  need  to  store 
the  forward  solution  in  order  to  form  the  adjoint  problem  and  the  cost  of  the  adjoint 
solve.  Our  approach  is  to  perform  the  adjoint  solves  using  relatively  coarse  scale 
discretizations  and  using  a  coarse  scale  representation  of  the  forward  solution  to  form 
the  adjoint  problem,  which  reduces  the  memory  overhead  and  the  cost  of  the  adjoint 
solve.  This  approach  is  motivated  by  the  following  observations. 

1.  Adjoint  problems  are  linear  and  often  present  fewer  numerical  difficulties  than 
the  associated  forward  problems. 

2.  Solutions  of  adjoint  problems  tend  to  vary  slowly  on  the  scale  of  the  dis¬ 
cretization,  whereas  residuals  of  forward  solutions  tend  to  oscillate  on  the 
scale  of  the  discretization 

3.  The  accuracy  required  of  the  adjoint  solution,  which  is  being  used  only  for 
error  estimation,  is  orders  of  magnitude  less  than  generally  desired  for  the 
forward  solution. 

An  enormous  literature  on  adaptive  methods  for  differential  equations  has  devel- 
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oped  over  nearly  six  decades  of  activity  and  the  major  developments  form  a  highly 
inter-connected  web.  We  do  not  attempt  to  review  the  history  of  adaptive  methods  or 
to  present  a  comprehensive  list  of  references.  Instead,  we  provide  only  a  short  list  of 
references  that  cither  contain  further  references  and/or  address  computational  issues 
related  to  adaptive  mesh  refinement  for  evolution  problems  [8,  7,  5,  4,  18,  22,  9,  10, 


17,  19,  15,  3,  1,  23,  24,  20,  2,  14]. 


This  paper  considers  adaptive  mesh  refinement  from  a  different  point  of  view  than 
much  of  the  existing  literature.  Namely,  we  are  concerned  with  trying  to  understand 
how  to  adapt  discretizations  based  on  under-resolved  solutions  on  relatively  coarse 
discretizations  in  order  to  obtain  particular  information,  as  opposed  to  analyzing 
adaptive  mesh  algorithms  in  the  asymptotic  limit  of  mesh  refinement.  This  point  of 
view  is  important  for  many  large  scale  applications,  for  which  such  conditions  are 
generic.  In  §2  we  review  the  standard  a  posteriori  error  imalysis  and  modify  this  for 
a  block  adaptive  strategy.  We  review  adaptive  error  control  in  §3  and  introduce  new 
features  necessary  for  block  adaptivity  and  several  block  adaptive  strategies.  One- 
and  three-dimensional  illustrative  computational  examples  are  provided  in  §4  and  we 
draw  conclusions  in  §5. 

2.  Discretization  and  error  estimation.  We  begin  by  reviewing  discretiza¬ 
tion  and  a  posteriori  error  estimation  for  evolution  problems  and  then  describe  the 
block-wise  discretization  and  present  the  corresponding  error  estimate. 

2.1.  Discretization.  We  formulate  the  discretization  as  a  space-time  finite  ele¬ 
ment  method  because  that  is  convenient  for  deriving  a  posteriori  error  estimates  based 
on  variational  analysis.  However,  we  emphasize  that  the  estimates  can  be  extended 
to  a  wide  range  of  discretizations,  e.g.  finite  difference  and  finite  volume  methods, 
which  can  be  written  as  equivalent  finite  element  methods. 

We  describe  two  finite  element  space-time  discretizations  of  (1.1)  called  the  con¬ 
tinuous  and  discontinuous  Galerkin  methods,  see  [11,  13,  12,  10,  17,  15].  We  partition 
[0,  T]  as  0  =  fo  <  ti  <  <2  <  •  •  ■  <  tri  <  •  •  ■  <  tiv  =  T,  denoting  each  time  interval  by 
In  =  (tri-iitn]  Aiul  time  Step  by  A'„  =  -  tn-i  and  we  construct  a  discretization  T 

of  fl  such  that  the  union  of  the  elements  in  7”  is  fl  while  the  intersection  of  any  two 
elements  is  either  a  common  edge,  node,  or  is  empty.  We  assume  that  the  smallest 
angle  of  any  element  is  bounded  below  by  a  fixed  constant.  To  mei^suro  t  he  size  of  the 
elements  of  T,  we  use  a  piecewise  constant  function  h,  the  so-called  mesh  function, 
defined  so  =  diam(A)  for  A  6  T.  Similarly,  we  use  k  to  denote  the  piecewise 
constant  function  that  is  on 

'I  he  approximations  are  polynomials  in  time  and  piecewise  polynomials  in  space 
on  each  space-time  “slab”  =  fl  x  In  space,  we  let  V  C  denote  the  space 

of  piecewise  linear  continuous  functions  defined  on  T,  whore  each  function  is  zero  on 
DO..  Then  on  each  slab,  we  define 


Finally,  we  let  denote  the  space  of  functions  defined  on  the  space-time  domain 
$2  X  [0,T]  such  that  u]^^  €  for  n  >  I.  Note  that  functions  in  W’’  may  be 
discontinuous  across  the  discrete  time  levels  and  we  denote  the  jump  across  by 
[te]„  =  -  w-  where  =  liin,,_,t„t  u;(s). 

We  use  a  projection  operator  into  F,  Pv  6  V,  e.g.  the  iJ  projection  satisfying 
{Pv,w)  =  {v,w)  for  all  w  C  V,  where  (■,•)  denotes  the  £2(^2)  inner  product.  We 
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use  the  II  II  for  the  La  norm.  We  also  use  a  projection  operator  into  the  piecewise 
polynomial  functions  in  time,  denoted  by  7r„  :  L^(/„)  where  P’(/n)  is  the 

space  of  polynomials  of  degree  q  or  less  rlofinod  on  /„.  The  global  projection  operator 
TT  is  defined  by  setting  tt  =  7r„  on  Sn- 

Definition  2.1.  The.  discontinuous  Galerkin  dG(q)  approximation  U  C  W’' 
satisfies  C/J"  =  Puq  and 


for  all  V  e  W®,  1  <  n  <  W  (2.1) 


We  also  use  a  related  method  for  solving  the  adjoint  problem: 

Definition  2.2.  The  continuous  Galerkin  cG(q)  approximation  U  £  W’ 
satisfies  Uq  =  Puq  and 


/ 


r  + (eVf/,Vu))d<=  (/(f/),r)d< 

1  ''tu-  I 

for  all  t>  £  I  <  n  <  N, 


(2.2) 


Note  that  U  is  continuous  across  time  nodes  when  the  space  mesh  is  fixed. 

With  appropriate  use  of  quadrature  to  evaluate  the  integrals  in  the  variational 
formulation,  those  Galerkin  methods  yield  standard  difference  schemes.  If  1,he  lumped 
mass  quadrature  is  used  in  space,  then  the  discrete  system  yielding  the  dG(0)  approxi¬ 
mation  is  the  same  as  the  system  obtained  for  the  nbdal  values  of  the  “backward  Euler 
in  time”- “second  order  centered  differenee  scheme  in  space”  finite  difference  scheme. 
Likewise,  the  cG(l)  method  is  related  to  the  Crank-Nicolson  scheme,  and  the  dG(l) 
method  is  rclatoil  to  the  third  order  sub-diagonal  Fade  difference  scheme.  Under  gen¬ 
eral  assumptions,  the  cG(q)  and  dG(q)  have  order  of  accuracy  g  +  1  in  time  at  any 
point  and  a  superconvergcnco  order  of  2g  -F  1  and  2q  respectively  at  time  nodes. 

2.2.  An  o  posteriori  error  estimate.  We  begin  by  defining  a  suitable  adjoint 
problem  for  error  analysis.  A  more  detailed  description  is  given  in  [15].  The  adjoint 
problem  is  a  parabolic  problem  with  coefficients  obtained  by  linearization  around  an 
average  of  the  true  and  approximate  solutions. 


(2.3) 


The  regularity  of  u  and  U  typically  imply  that  /  is  piecewise  continuous  with  respect 
to  t  and  a  continuous,  function  in  spac:e. 

Written  out  pointwise  for  convenience,  the  adjoint  problem  to  (1.1)  for  the  gen¬ 
eralized  Green’s  function  associated  to  the  data  rp,  which  determines  the  quantity  of 
interest. 
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is 

-  V  ■  {(V(p)  -  J(j)  -  xj),  (x,t)  G  ft  X  (T,0], 

<  4)(x,  t)  =  0,  (x,  t)  fe  an  X  (7\  0],  (2.4) 

((^(x,7’)  =  0,  X  €  n, 

This  choice  for  the  adjoint  yields  the  following  error  representation  formula  for 
the  dG  method. 

Theorem  2.3.  dG  A  Posteriori  Error  Estimate 

/  ie,rP)dt  =  ((/  -  P)uo,m)  T  ^{[UU.u{^P4>~<P)t-i) 

v  0  _ 1 

nssi 

+  f  {{U,7rP(l>-d>)  +  {t{U)VU,WinP<i>-^))  -{f{U),TrP(P-(P))dt.  (2.5) 

J  0 

The  initial  error  is  e“(0)  =  (/  —  P)uo. 

In  practice,  we  compute  a  numerical  solution  of  the  linear  adjoint  problem  ob¬ 
tained  from  (2.4)  by  replacing  a  with  the  computed  approximate  solution  U  in  the 
definition  of  /  and  solve  using  a  higher  order  method  in  space  and  time,  see  [15].  We 
denote  the  approximate  adjoint  solution  by  <E>.  We  focus  on  the  dG  method,  while 
application  to  the  cG  method  is  analogous. 

Corollary  2.4.  The  approximate  a  posteriori  error  estimate  for  the  dG  method 
is 


\fi. 


xp)  dt 


E{U)=^E{U-,xp) 


N 


((/-P)uo,$(0))  +  ^([t/]„_,,(,rP$-$)+_i) 


n=l 


^  {{U,  ttP^  -  4>)  +  (£([/) Vt/,  V(7rP$  -  $))  -  (f{U),  7rP4>  -  $))  dt 


(2.6) 


2.3.  Blockwise  discretization.  We  describe  the  blockwise  formulation  of  the 
discontinuous  Galerkin  method.  We  partition  [0,  Tj  into  time  blocks  0  =  7o  <  Ti  < 
T2  <  ■■■  <Tb  <  ■•■  <Tb  —  T.  We  discretize  each  block  [Tb-i,Th]  by  Tb-i  =  tb,o  < 
<6,1  <  ■  •  •  <  tb,Nh  =  Tb,  denoting  each  subinterval  by  =  (tb.n-i,  4,n]  and  time  step 
f^h.n  =  tb,n  —  tb,n-i-  To  cach  block  [7fc_i,r(,],  we  associate  a  discretization  %  of 
it  arranged  so  the  union  of  the  elements  in  Th  is  while  the  intersection  of  any  two 
elements  is  either  a  common  edge,  node,  or  is  empty.  We  assume  that  the  smallest 
angle  of  any  element  is  bounded  below  by  a  fixed  constant.  To  measure  the  size  of 
the  elements  of  %,  we  use  the  mesh  function  hb. 

The  approximations  are  polynomials  in  time  and  piecewise  polynomials  in  space 
on  each  space-time  ^‘slab”  5fc,„  =  fi  x  In  space,  we  let  14  c  denote  the 

space  of  piecewise  linear  continuous  functions  defined  on  %,  whore  each  function  is 
zero  on  dQ.  Then  on  each  slab,  we  define 

XM{x,t)  :  xv{x,t)  =  Y^PVb.jix),  VbJ  e  14,  (X,t)  e  Sb,n 
j=o 

Finally,  we  let  W  denote  the  space  of  functions  defined  on  the  space-time  domain 
it  X  [0,7’]  such  that  u|si,,„  €  for  6,  n  >  1.  Note  that  functions  in  may  be 
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discontinuous  across  the  discrete  time  levels  and  we  denote  the  jump  across  by 
[u'lfc.n  = 

To  compute  the  dG  approximation  on  the  new  block,  we  project  the  final  value  of 
the  approximation  from  the  previous  block  onto  the  new  mesh.  We  use  a  projection 
operator  Pi,v  £  Vj,  and  a  projection  operator  into  the  piecewise  polynomial  functions 
in  time,  denoted  by  ->  We  then  define  ni,  as  ni,  =  7rt,.„  on 

Sb,n-  Finally,  we  define  global  projections  P  and  tt  which  on  each  block  are  Pb  and 
rcb  respectively. 

Definition  2.5.  The  blockxvise  discontinuous  Galerkin  dG(q)  approximation  U  £ 
W’’  satisfies  =  P\uq  and  for  6  =  1, 2,  •  ■  •  ,B, 

/  {{U,v)  +  {fVU,Vv))dt+ {[U]b,„-i,v^)  -  [  {f{U),v)dt 

for  all  V  €  ,  1  <  n  <  iVfc.  (2.7) 


2.4.  A  blockwise  a  posteriori  error  estimate.  Adapting  the  standard  argu¬ 
ment  that  yields  (2.5),  we  obtain  a  blockwise  a  posteriori  error  e.stimate. 

Theorem  2.6.  Blockwise  A  Posteriori  Error  Estimate 


T  a 

{e,rlr)dt  «  ((/  -  Po)no,$(0))  +  ^((/  -  Pb)U,^Tb-i)) 

h=l 

^  /  rn 

+  Z,  (  /  +  (Kt^)V(7,  V(^n$  -  <!>))  -  (/(t/),  rtPb^  -  $))  dt 

6=1  '  ■'^1-1 

+  Z  ([t^kn-l,(7rm-$)6Vl 

n-  =  l 

The  second  term  on  t.hc  right  measures  the  effects  of  changing  meshes  on  the  accuracy 
of  the  approximation.  A  similar  “jump”  term  already  appears  in  the  estimate  for  the 
standard  dG  method  at  each  time  step.  In  this  case  of  transitions  between  blocks,  the 
“jump”  arises  because  of  mesh  changes  between  blocks.  Note  that  the  adjoint  weight 
does  not  involve  the  projection  of  $  into  the  approximation  space  (i.e.  Galerkin 
orthogonality).  Instead,  the  contributions  from  the  projections  accumulate  in  the 
same  way  as  an  initial  error. 

Our  purpose  is  to  use  the  a  posteriori  bounds  to  choose  block  times  {Tb} 

and  corresponding  meshes  %  and  timesteps  fcb.i-  An  important  issue  is  the  effect  of 
transferring  solutions  between  the  meshes  of  adjacent  blocks  on  the  accuracy  of  the 
computed  information,  and  so  we  address  the  computation  of  a  bound  on  the  second 
term  on  the  right  in  (2.8), 


a 

(2.9) 

6=1 

3.  Adaptive  error  control.  We  start  off.  by  describing  some  standard  ap¬ 
proaches  to  adaptive  error  control  and  the  relation.to  adaptive  error  control  based  on  a 
posteriori  error  estimates.  We  then  turn  to  the  problem  ‘of  choosing  blocks  for  a  block 
discretization  and  generating  the  corresponding  spatial  and  temporal  discretizations 
for  each  block. 
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3.1.  Goal  oriented  adaptive  error  control.  The  aim  of  goal  oriented  adap¬ 
tive  error  control  Ls  to  generate  a  mesh  with  a  nearly  minimal  number  of  elements 
such  that  for  a  given  tolerance  TOL  and  data  Vh 


(3.1) 


We  note  that  (3.1)  cannot  be  verified  in  practice  because  the  error  is  unknown,  so 
instead  we  use  an  estimate  or  a  bound  for  the  error  in  the  quantity  of  interest.  Different 
ways  to  generate  an  acceptable  mesh  vary  by  the  estimate  or  bound  used  for  the 
(}uantity  of  interest  as  well  as  the  strategy  for  mesh  refinement. 

For  example  using  the  a  posteriori  estimate  (2.6),  the  goal  of  adaptive  error  control 
is  to  determine  a  discretization  so  that  a  mesh  acceptance  criterion. 


EiU)  <  TOL,  (3.2) 

is  satisfied.  If  (3.2)  is  not  satisfied,  then  we  refine  the  mesh  in  order  to  compute  a  new 
solution  for  which  the  criterion  is  met.  Refinement  decisions  require  identifying  the 
contributions  to  the  error  from  discretization  on  each  element.  We  can  write  E{U)  as 
a  sum  over  space-time  elements, 


E{U) 


((/  -  ^)«o,^(o))a  +  e  e 


A€T 


=1  Aer 


+  E  E  /  ((f>.’rP$-$)A  +  (e(f/)VI7,V(7rP$-$))^-(/(r/),7rP$-«>)A)df 

n=i  Aer 


where  (  ,  )a  denotes  the  inner  product  on  element  A.  This  clearly  identifies 
possible  element  contributions.  , 

However,  a  major  difficulty  is  that  the  error  estimate  generally  involves  a  large 
amount  of  cancelation  among  the  element  contributions,  which  makes  determining  a 
truly  effic:icnt  refinement  strategy  extremely  difficult. 


Fig.  3.1.  The  element  contributions  to  the  error  in  integration. 


Example  3.1.  We  consider  a  first  order  accurate  numerical  solution  that  has  the 
element  contributions  shown  in  Fig.  3.1. 

The  first  time  step  has  the  largest  contribution.  The  next  three  steps  each  con¬ 
tribute  -0.033,  so  cancelation  means  that  the  total  contribution  from  the  first  four 
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steps  is  0.001.  Likewise,  the  next  six  steps  contribute  +0.003  in  total.  The  last  four 
steps  contribute  0.08  in  total.  The  total  error  is  therefore 

.1  -  3  X  .033  +  .011  -  .01  +  .011  -  .01  +  .011  -  .01  +  4  X  .02  =  0.084 

If  we  use  a  standard  approach  of  refining  only  some  fraction  of  the  elements  with  the 
largest  contributions,  we  are  likely  to  refine  the  first  four  steps.  For  simplicity,  we 
assume  that  the  elements  marked  for  refinement  arc  divided  into  two  time  steps.  The 
resulting  integration  will  have  accuracy 

^x2x.l-ix6x  .033  f  .011  -  .01  f  .011  -  .01  +  .011  -  .01  +  4  x  .02  w  0.0835. 

Note  that  the  individual  clement  contributions  decrease  at  a  second  order  rate.  The 
problem  is  that  even  though  the  element  contributions  in  the  first  four  steps  are 
individually  large,  there  is  significant  cancelation  and  refinement  in  this  region  and 
refinement  does  not  decrease  the  error  significantly.  On  the  other  hand,  if  we  refine 
the  last  four  time  steps  instead,  we  obtain 

.1  -  3  X  .033  +  .011  -  .01  +  .011  -  .01  +  .011  -  .01  +  ^  X  8  X  .02  w  0.044. 

While  this  is  a  non-standard  approach,  it  decreases  the  error  significantly. 

In  the  adjoint-weight  approach,  the  issue  of  cancelation  of  error  is  neglected  in  a 
sense  by  replacing  the  accurate  error  estimate  E{U)  by  an  inaccurate  upper  bound, 

E{U)<M{U)^M{U-,^),  (3.3) 

where  we  define  iE  ({/;  ^)  by  summing  bounds  over  each  element. 

Definitio.n  3.2.  Element-wise  upper  bound  on  the  total  error 


N 


^{u-,rp)=  + 

ACT  n=lAer 

n=lA6T  •'E.-i 


dt 


Thus,  if  (3.2)  is  not  satisfied,  the  mesh  is  refined  in  order  to  achieve 

/E(f/)<TOL.  (3.4) 

The  adaptive  error  control  problem  can  now  be  profitably  posed  as  a  constrained 
minimization  problem,  namely  to  find  a  mesh  with  a  minimal  number  of  degrees  of 
freedom  on  which  the  approximation  satisfies  the  bound  (3.4).  Using  the  fact  that 
the  bound  jE  is  a  sum  of  positive  terms  and  assuming  tlu;  .solution  is  asymptotically 
accurate,  a  calculus  of  variations  argument  yields  the  generic  (see  e.g.  [9,  10,  3,  2]). 

Principle  of  Equidistribution  An  approximate  solution  of  the 
constrained  optimization  problem  for  an  optimal  mesh  for  an  upper 
bound  on  the  error  is  achieved  when  the  elements  contributions  to 
the  bound  are  approximately  equal. 
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The  Principle  of  Equiciistribution  has  been  used  in  various  forms  at  least  since  the 
seventies  (and  probably  earlier  in  industry).  However,  experience  with  a  wide  range 
of  probhuns  suggest  that  the  bound  IE  {U)  is  generically  several  orders  of  magnitude 
larger  than  the  estimate  E{U).  A  strategy  based  on  the  Principle  of  Equidistribution 
that  optimizes  computational  cost  with  respect  to  a  error  bound  and  not  the  actual 
error  can  therefore  result  in  significant  over- refinement. 

In  general,  there  are  many  solutions  of  the  constrained  minimization  problem 
associated  with  (3.4).  An  adaptive  mesh  algorithm  is  a  procedure  for  computing  an 
acceptable  solution.  Traditionally,  different  approaches  are  used  for  spatial  and  tem¬ 
poral  adaption.  A  global  “compute-estimate-mark-adapt”  algorithm  (see  for  example 
1.1)  is  typically  used  for  spatial  meshes.  This  is  an  iterative  approach  in  which  only 
some  fraction  of  the  elements  on  which  the  contribution  to  the  error  bound  is  largest 
are  refined  during  each  iteration  and  whole  cycle  is  iterated  until  a  prescribed  tol¬ 
erance  is  achieved.  Temporal  approaches  to  mesh  adaption,  e.g.,  local  error  control 
[21],  tend  to  use  a  “sweeping”  strategy  from  initial  to  final  time,  where  a  solution 
is  advanced  past  each  time  step  only  when  the  step  contribution  is  estimated  to  be 
lower  than  an  acceptable  fraction  of  the  total  error.  This  may  be  viewed  sis  a  gener¬ 
ally  pessimistic  way  to  achieve  the  Principle  of  Equidistribution  because  it  removes 
positive  effc’cts  of  cancelation  of  ('rror  altogether.  As  a  consequence  of  those  differ¬ 
ences,  element  contributions  to  the  error  estimate  or  bound  typically  vary  in  size  quite 
considerably  while  contributions  from  different  time  intervals  arc  more  nearly  equal. 

We  use  a  strategy  that  treats  space  and  time  discretizations  more  equitably.  In 
the  case  of  a  parabolic  problem,  it  is  straightforward  to  distinguish  the  time  and  space 
contributions  to  the  bound  IE .  We  define  the  time  and  space  bounds  as  follows. 

Definition  3.3.  Element-wise  temporal  and  spatial  error  bounds 


N 


n=i  Aer ' 


+  E  E  /  +  (£(t/)Vt/,  v(7r  -  /)P$)a 

n=i  Aer  ''*'•  -1 


-{f{U),{n-I)P^)^dt 


,  (3.5) 


^z{U)=  |((^-^H,$(0))a|  f-E  E  ”  ([^^1n-i.(P$-^)+_i) 

■T 

N  I 

T.Z  .  {U,  P4)  -  0)a  f  (f  ([/)  Vt/,  V(P'I>  -  4>))a 

ri=lA6r'‘^‘"  I 


Aer 


+ 


A 


-(/(f/),P$-$)Adt 


.  (3.6) 


We  sec  that  the  time  bound  is  precisely  the  a  posteriori  bound  for  the  dG  approxi¬ 
mation  for  the  “method  of  lines”  initial  value  problem  resulting  after  discretization  in 
space.  The  adjoint  weight  depends  on  the  projection  of  the  adjoint  solution  into  the 
time  finite  element  space.  On  the  other  hand,  the  adjoint  weight  in  the  space  bound 
depends  on  the  projection  of  the  adjoint  solution  into  the  spatial  finite  element  space. 
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Wo  split  the  error  between  the  time  and  space  contributions  and  refine  the  current 
mesh  in  order  to  achieve 

'rni  TOT 

^  —  and  £,{U)  <  ^  .  (3.7) 

On  a  given  time  interval,  this  requires  an  iteration  during  which  both  the  space  mesh 
and  time  steps  are  refined. 

3.2.  Goal  oriented  block  adaptive  error  control.  For  the  purpose  of  de¬ 
veloping  a  block  adaptive  algorithm,  we  treat  adaptivity  with  respect  to  space  and 
time  in  the  same  way.  The  reason  is  that  we  determine  the  blocks  by  predicting  the 
local  element  sizes  (or  number  of  sub-elements)  that  are  required  in  the  final  mesh. 
We  create  a  block  by  grouping  together  a  set  of  coarse-scale  space-time  slabs  that  are 
adjacent  in  time  and  satisfy  some  criteria,  e.g.  similar  spatial  meshes  are  predicted 
for  the  space-time  slabs  in  the  block  or  a  maximal  number  of  elements  are  predicted 
to  be  required  in  the  block. 

3.2.1.  Choosing  a  global  tolerance  for  the  error  bound.  We  want  the  pre¬ 
dictions  of  the  element  sizes  required  in  an  acceptable  fine  scale  mesh  to  be  as  accurate 
as  possible.  We  recall  that  an  acceptable  mesh  need  only  satisfy  the  estimate  criterion 
(3.2)  and  not  the  more  stringent  bound  criterion  (3.4).  We  define  the  overcstimation 
factor  for  a  given  mesh. 


’  E{U)  ’ 

and  the  corres])onding  absolute  tolerance  for  , 


ATOL  =  7  X  TOL . 


We  replace  (3.4)  by 


(3.8) 


Note  that  ATOL  ss  TOL  when  there  is  little  cancelation  among  the  element  con¬ 
tributions  and  ATOL  >  TOL  otherwise.  In  this  way,  we  attempt  to  mitigate  the 
inefficiency  that  is  introduced  by  replacing  an  accurate  error  estimate  by  an  inaccu¬ 
rate  bound  in  decisions  about  mesh  refinement.  This  approach  for  setting  tolerances 
is  discussed  further  in  [16]. 

3.2.2.  Predicting  refinement  in  space.  Given  a  local  space-time  element  6  = 
fc5(A,ri)  =  A  X  in  the  spacc'-time  slab  that  is  marked  for  refinement,  we 

show  how  to  predict  the  number  of  space-time  elements  that  are  needed  to  meet  the 
acceptance  criterion.  We  assume  that  in  the  current  mesh,  there  are  N  time  steps  and 
M  space  elements  in  each  space-time  slab,  giving  a  total  of  NM  space-time  elements. 
We  define  a  local  absolute  tolerance 


LATOL  = 


ATOL 
2NM  ■ 


By  the  Principle  of  Equidi.stribution,  we  adopt  the  goal  of  refining  each  space-time 
element  so  that  the  local  element  contribution  is  approximately  LATOL. 
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Using  a  priori  convergence  analysis,  see  [15],  it  is  possible  to  show  that  there  is  a 
constant  C  such  that 


(3.9) 

as  0,  where  p  is  reflated  to  the  order  of  the  finite  element  met  hod  in  space  and 

/lA  is  the  element  size.  Likewise,  we  can  show  constant  C  such  that 

(3-10) 

as  A:  0,  where  q  is  related  to  the  order  of  the  finite  element  method  in  time. 

Now  suppose  that  an  element  6,h.w  in  the  final  mesh  is  obtained  from  ©oid  in  the 
current  mesh  by  refinement.  We  have 


LATOL  %  M 


E. 


'  Sold 


(3.11) 


This  yields  a  prediction  for  the  new  mesh  size. 


,  ( LATOL  \ 

^A„ow  ~  1  I  j  ^  ^Aoid-  (3.1.2) 

Recalling  that  d  is  the  space  dimension,  this  predicts  that  the  element  Aou  should 
be  refined  into  roughly 


sul)-elcments. 

3.2.3.  Predicting  refinement  in  time.  For  refinement  in  time, 

yE(|  wjEtL  X  «LATOL. 

IS„.,„  I  Sold  koM  J 

This  yields  a  prediction  for  the  new  mesh  size, 


/ latol\ 


1/9 


X  k 


old* 


This  predicts  that  the  time  step  /Coid  should  be  refined  into  roughly 


_  f  I  Sold 


1/9 


I LATOL 


(3.13) 


(3.14) 


(3.15) 


(3.16) 


sul>intervals. 
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3.2.4.  Determining  overlaps  for  meshes  on  adjacent  blocks.  After  the 
meshes  for  each  block  are  determined  based  on  the  a  posteriori  prediction  of  error,  we 
need  to  estimate  the  effects  of  transferring  the  solution  between  meshes  on  adjacent 
blocks.  See  §  4.1  for  an  example  that  illustrates  this  point.  Recall  that  (2.9)  provides 
a  bound  on  these  effects.  I'hc  difficulty  with  using  (2.9)  is  that  we  do  not  have  the 
fine  scale  numerical  solution  U  required  for  that  expression  until  after  solving  on  the 
fine  scale,  whereas  ideally  wo  could  predict  a  reasonable  overlap  before  (imputing  the 
expensive  fine  scale  solution. 

We  list  three  strategies  for  mitigating  the  possibility  of  projection  error  in  our 
block  adaptive  framework. 

1.  There  is  a  very  simple  strategy.  In  forming  the  space  mesh  for  the  block 
\Tb-i,Th]  X  12,  we  guide  refinement  by  using  the  maximum  of  the  element 
contributions  on  each  individual  element,  taking  the  mfvximum  over  the  time 
intervals  included  in  the  block.  We  may  simply  include  the  maximum  over 
the  last  time  interval  included  in  the  previous  block,  [V),.  2, i.e.,  over 
the  interval  [<6-1, Nt  i-i,<6-i,N(,  J-  We  can  be  even  more  conservative  by 
including  some  number  of  the  last  time  steps  in  the  maximum  computation. 

2.  We  can  use  gradient  recovery  [6]  to  compute  an  approximate  solution  on  the 
fine  scale  mesh  in  each  block  using  the  solution  from  the  last  time  interval 
contained  in  each  block.  We  can  then  directly  compute  (/  —  Fb)U  for  each  b 
and  evaluate  (2.9). 

3.  We  can  evaluate  (2.9)  a  posteriori  by  evaluating  {I  —  Ph)U  using  the  fine  scale 
forward  solution  and  the  coarse  scale  adjoint  solution. 

3.3.  Block  adaptive  algorithms.  Using  the  development  above,  we  present  a 
generic  block  adaptive  algorithm  in  Algorithm  3.1.  We  provide  a  detailed  algorithm 
in  Appendix  A. 

Algorithm  3.1  Block  Adaptive  Algorithm 

1:  Choose  the  “coarse”  mesh  and  time  step 

2:  Compute  the  coarse  .scale  numerical  solution 

3;  Estimate  the  element  contributions  to  the  error  for  the  current  solution 

4:  Predict  the  number  of  space-time  elements  into  which  each  current  space-time 
element  is  to  be  divided  using  (3.13)  and  (3.16) 

5:  Build  block  discretizations  by  constructing  meshes  satisfying  the  requirements  for 
groups  of  neighboring  time  steps 

6:  Compute  the  fine  scale  numerical  solution  using  the  block  discretizations 


We  note  that  the  Block  Adaptive  Algorithm  3.1  can  be  iterated,  so  that  the  fine 
scale  becomes  the  new  coarse  scale,  and  a  new  fine  scale  is  subsequently  computed. 
In  crude  terms,  the  block  adaptive  Algorithm  3.1  is  analogous  to  the  core  estimate- 
inark-refine  algorithm  at  the  heart  of  the  generic  Algorithm  1.1,  but  is  different  in 
the  mark  and  refine  steps.  The  critical  step  defining  the  block  adaptive  algorithm 
Algorithm  3.1  is  the  strategy  used  to  create  block  discretizations.  Once  the  blocks 
are  identified,  we  can  use  any  adaptive  mesh  refinement  strategy  for  producing  the 
actual  meshes.  We  describe  several  strategies  for  determining  block  discretizations. 

3.3.1.  A  memory-bound  strategy.  In  the  first  strategy,  we  assume  there  is 
a  target  number  of  elements  Nmax  in  space  that  is  maximal  in  some  sense,  e.g.,  the 
largest  number  of  elements  that  can  be  stored  in  core.  We  form  blocks  by  creating  a 
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union  of  adjacent  coarso-scale  space-time  slabs,  one  slab  at  a  time,  until  the  projected 
space  mesh  for  the  block  uses  Nmax  elements.  To  create  the  block  mesh,  we  use  the 
maximum  of  the  predicte<l  number  of  elements  Nelem.children  on  each  individual 
element  (given  by  equation  (3.13))  in  the  union  forming  the  block.  We  illustrate  in 
Fig.  3.2.  The  parameter  6  governs  how  often  the  mesh  is  replaced  by  a  coarser  mesh, 
where  ^  w  10  works  well  in  practice. 


Division  into  Space-Time  Blocks 


u 

S 

P 


Space 


Predicted  Mesh  for  the  First  Block 


Space 


Fia.  :i.2.  The  memory  bound  strategy  is  used  for  a  traveling  pulse  that  moves  with  constant 
speed  from  left  to  right.  Left:  The  original  uniform  mesh  and  a  contour  plot  of  the  mimber  of 
predicted  elements  of  new  sub-elements  Nelem.children.  The  scale  is  from  dark  (low)  to  white 
(high).  Right:  The  predicted  number  of  new  sub-elements  Nelem-children  for  the  first  block,  which 
consists  of  three  adjacent  space-time  slabs  from  the  original  discretization. 


3.3.2.  A  correlation  strategy.  In  the  second  strategy,  we  aim  to  choose  blocks 
in  order  to  use  a  relatively  small  number  of  elements,  so  Nmax  may  be  considerably 
smaller  than  for  the  first  algorithm.  This  strategy  forms  a  block  by  grouping  to¬ 
gether  adjacent  coarse-scale  space-time  slabs  whose  predicted  number  of  elements 
Nelem_children  are  close. 

In  [14],  we  consider  the  problem  of  detecting  significant  overlap  of  local  element 
contributions  for  different  computations.  Following  the  approach  there,  given  two 
vectors  v,w  whose  coefficients  are  element  contributions  to  an  error  estimate,  we 
define  their  correlation  to  be  c{v,  w)  =  v  ■  w.  We  say  that  v  is  significantly  correlated 
with  w  if 


c{v,  w) 
||ui||2 


>  7i  and 


II  a; 


<  72- 


where  0  <  71,72.  The  first  condition  insures  that  v  has  a  suitable  large  projection 
onto  vi  while  tlu^  second  condition  corrects  for  differences  in  scale  between  v  and  w 
(con-sider  ||i7||  »  ||ii;||  so  that  c(v,w)  »  ||ii)||). 

We  implement  the  new  criterion  for  creating  blocks  by  choosing  to  add  the  next 
time  slab  to  a  current  block  based  on  the  correlation  criterion. 

3.3.3.  Global  strategies.  In  the  first  two  strategies  for  creating  blocks,  we 
sweep  through  time.  We  can  also  use  a  bisection  search  beginning  with  the  original 
large  block  and  subdividing  to  find  acceptable  blocks.  In  analog  to  the  difference 
between  the  standard  global  strategy  for  space  mesh  refinement  to  achieve  the  Prin¬ 
ciple  of  Equidistribution  and  the  local-error  control  approach,  the  bisection  search  is 
a  global  strategy  that  can  be  a  more  efficient  way  to  achieve  equidistribution. 
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4.  Computational  Examples.  Wc  apply  the  block  adaptive  algorithms  to  sev¬ 
eral  prototypical  examples  in  one  and  three  space  dimensions.  The  one  dimensional  ex¬ 
amples  illustrate  several  key  points  when  implementing  block-adaptive  methods,  while 
the  three  dimensional  examples  include  a  traveling  wave  front,  a  solution  that  under¬ 
goes  time-  and  space-localized  perturbations,  and  a  periodic  motion  in  a  convection- 
dominated  flow.  i 

The  forward  problems  and  adjoint  problems  arc  solved  with  linear  and  quadratic 
elements  in  space  and  dGO  and  cGl  in  time  respectively.  The  one  dimensional  ex¬ 
amples  are  computed  using  the  Matlab  code  ACES  [25].  The  three  dimensional  ex¬ 
amples  are  performed  on  a  hexaliedral  mesh  using  a  trilincar  spatial  basis  for  the 
forward  problem  and  a  triquadratic  basis  for  the  adjoint.  Local  mesh  refinement  is 
accomplished  by  the  use  of  hanging  nodes  where  one  hanging  node  per  edge  or  face 
is  allowed.  Conformity  of  the  basis  is  obtained  by  interpolation  of  the  surromiding 
regular  nodes.  Ihe  use  of  an  hierarchical  octree-based  data  structure  assists  refine¬ 
ment  but  also  allows  for  de-refinement  when  the  element  indicators  are  small.  For  the 
convection  driven  flow  problem,  SUPG  is  employed  for  both  the  forward  and  adjoint 
problems,  with  piu-ameter 


(1/Af  +  Ujh)  ' 

where  At  is  the  time  step  and  U  is  the  speed  of  the  convection  field  at  the  tmrrent 
time,  i.e.,  U  =  ||/3||2  in  (4.5).  This  is  not  an  obstacle  for  the  block-adaptive  frame¬ 
work,  as  we  simply  modify  the  theoretical  convergence  rate  p  in  the  computation  of 
Nelem_children  in  (3.13). 

4.1.  Example  One:  Projection  errors  between  blocks.  VVe  illustrate  the 
necessity  for  addressing  the  effect  of  transferring  solutions  between  space-time  blocks 
with  a  simple  one-dimensional  example  involving  a  traveling  wave. 

ui  -  Uxx  =  /(a:,t),  0  C  X  <  1,  0  <  (, 

<  u(0,t)  =  u(l,f)  = /3(f),  0<f,  (4.1) 

(i(.r,())  =  tanh(o:(x  -  0.2)),  0  <  x  <  1, 

where  a  =  50  and  /  and  ft  are  chosen  to  give  an  exact  solution  u  =  tanh(a(x-f-0.2)). 
We  solve  with  a  coarse  mesh  using  /i  =  0.1  and  time  step  k  =  0.05  from  initial  time 
0  to  final  time  0.6.  The  quantity  of  interest  is  the  average  space-time  error.  We 
compute  a  fine  scale  solution  using  two  blocks  derived  from  the  coarse  scale  solution. 
The  first  block,  t  =  [0,0.3],  uses  a  finer  spatial  mesh  in  the  region  x  e  [0.1, 0.6],  while 
the  second  block  uses  a  fine  mesh  in  the  region  [0.5, 1],  so  the  overlap  is  minimal  and 
and  the  prcdicticms  for  refinement,  areas  are  incorrect.  Consequently,  the  approximate 
traveling  wave  travels  too  quickly.  The  first  block  solution  at  t  =  0.3  and  its  projection 
onto  the  .second  block  at  t  —  0.3  is  displayed  in  Fig.  4.1. 

In  1*  ig.  4.1  we  illustrate  the  a  posteriori  use  of  (2.9)  to  correct  the  projection 
error.  Block  1  is  computed  using  the  predicted  fine  scale  mesh.  Block  2  is  tested  for 
significant  projection  error  using  (2.9)  using  the  fine  scale  solution  for  Block  1  and 
the  mesh  for  Block  2  is  refined  if  the  elementwise  projection  error  exceeds  LATOL. 
We  note  that  the  overlap  strategy  for  the  projection  error  in  §3.2.4  also  works  well  in 
this  particular  example. 

4.2.  Example  Two:  Coarse  scale  resolution.  Since  wc  are  using  the  coarse 
scale  discretization  to  predict  the  global  behavior  of  the  solution  on  the  fine  scale. 
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Fig.  4.1.  Problem  H-l).  The  circles  indicate  the  spatial  meshes  used  in  each  of  the  two  blocks. 
Left:  the  solution  on  Ulock  I.  Middle:  the  projection  of  the  approximate,  solution  in  block  1  onto 
the  mesh  in  block  2.  Right:  the  solution  onto  Block  2  after  using  the  projection  error  estimate 
(2.9)  to  correct  signifir/int  projection  errors  between  the  two  blocks.  This  demonstrates  the  possible 
consequences  when  the  meshes  for  neighboring  blocks  do  not  overlap  sufficiently. 


it  is  important  to  insure  that  the  coarse  scale  discretization  is  not  too  coarse.  (This 
is  a  difference  between  the  block  adaptive  approach  and  a  standard  adaptive  mesh 
refinement,  which  is  generally  started  with  a  very  coarse  mesh.)  This  issue  is  especially 
important  for  nonlinear  problems  since  linciarization  is  used  to  define  the  adjoint 
problem,  which  in  turn  provides  the  means  to  quantify  the  effects  of  cancelation  and 
accumulation  of  errors. 

Consider  the  one-dimensional  nonlinear  parabolic  equation 

iH  =  «(w  -  1)(1  n^),  -1  <  .r  <  1,  0  <  f  <  0.6, 

<  u(0,f)  = -1,  n(l,t)  =  1,  0<f,  (4.2) 

7/(.'r,0)  =  tanh(rv(.7;  -  0.2)),  -1  <  x  <  1, 

We  choose  a  to  obtain  the  same  solution  as  the  example  in  §  4.1,  u  =  tanh(a(x’  — 
t  —  0.2)).  The  quantity  of  interest  is  the  average  space-time  error.  For  the  coarse 
di.scretization,  we  use  h  =  0.05  and  k  =  0.05.  These  choices  provide  an  excellent 
coarse  scale  discretization  for  the  linear  example  in  §  4.1  but  docs  not  work  well  for 
the  nonlinear  version.  We  show  two  snapshots  of  the  solution  u  in  Fig.  4.2  at  t  =  0.3 
and  t  =  0.6.  The  wave-speed  is  predicted  inaccurately,  which  leads  to  a  poor  block 
selection  and  this  subsequently  affects  the  fine  scale  accuracy.  Using  a  (X)arse  scale 
discretization  with  h  =  0.1  and  A:  =  0.1  yields  inaccurate  results. 

I, 

0.5 
0 

-0.5 


1.5 1 - - - - - - - - 

0  02  0.4  0.6  0.8  I 

Fig.  4.2.  Problem  (4-2).  CoiTclatxon  strategy  with  an  insufficiently  accurate  coarse-scalc  solu¬ 
tion.  Solution  on  the  adapted  mesh  at  t  -  0.3  and  t  =  0.6  respectively. 


The  poor  predictions  based  on  the  coarse-scale  discretization  can  be  avoided  by 
slightly  enriching  the  discretization  with  a  finer  time  step.  We  use  a  coarse  discretiza- 
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tion  with  h  =  0.05  and  k  =  0.01  and  the  correlation  strategy  to  produce  blocks.  The 
approximate  solution  on  the  adapted  mesh  at  t  =  0.45  is  shown  in  Fig.  4.3. 


Fig.  Problem  (4-^)-  Correlation  strategy  unth  an  improved  coarse-scale  solution.  Solution 
on  the  adapted  mesh  at  t  =  0.45  on  blocks  3  and  4  Tespectiveiy. 


4.3.  Example  three:  A  traveling  w^ave  solution.  This  example  is  a  wave 
propagating  along  tlu'  main  diagonal  of  the  unit  cube  (0  =  [0, 1]  x  [0, 1]  x  [0, 1]).  The 
governing  equation  is 


Au  =  f{x,t.), 

<  u{x,t)  =  0, 

u{x,0)  =  (.7.1  -  7:‘f)(.7:2  -  .r2)(.7;3  -  -4)  arctan(^^3:f  +  .7;|  +  .7:|), 
where  c  =  75  and  /  is  constructed  to  yield  the  exact  solution 


X  c  il,  Q  <t, 

X  &  911, 0  <  t, 
X  e  ii, 

(4.3) 


VZ  .C\/3  /“7  y  ~ 
u  =  —  arctan(— —  yif  +  +  X3  -  t). 


:\/3 


The  coarse  block  solution  uc  is  constructed  on  an  8  x  8  x  8  uniform  mesh  using 
hexahedral  meshes  with  an  initial  time  step  of  0.1.  The  quantity  of  interest  is  the 
time  average  of  the  solution  value.  The  memory  bound  strategy  is  used  to  construct 
the  discretization  blocks  with  ATOL  =  0.000178  and  Nmax=50obo.  Block  information 
is  given  in  Table  4.1.  As  might  be  expected,  all  of  the  blocks  use  approximately  the 
same  number  of  elements.  We  show  contour  plots  of  the  solution  on  “slices”  of  some 
of  the  block  meshes  along  the  plane  x  =  0.5  in  Fig.  4.4. 


Block 

Tb-i 

Tb 

#  vertices 

#  hexahedra 

1 

0 

0.4 

58711 

50394 

2 

0.4 

0.6 

63219 

54503 

3 

0.6 

0.7 

72267 

61265 

4 

0.7 

0.8 

62626 

52368 

5 

0.8 

1 

64764 

54860 

6 

1 

1.1 

62790 

54377 

Table  4.1 

Problem  (4-3).  Blocks  resulting  from  the  memory  bound  strategy. 


4.4.  Example  Four:  Localized  forcing  in  space  and  time.  This  example 
contrasts  the  difference  in  the  blocks  produced  by  the  memory  bound  and  correlation 


blockwisp:  adaptivity 
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Fig.  4.4.  Problem  (4-3).  Memory  bound  strategy.  Slices  through  the  mesh  perpendicular  to  the 
x-axis  at  X  =  0.5.  Upper  left:  t  —  0  (block  }).  Upper  right:  t  =  0.44  (block  2).  Lower  left:  t  =  0.6 
(block  3).  Lower  right:  t  =  1.1  (block  6). 


strategies  when  solving  an  equation  with  source  terms  that  are  localized  in  space  and 
time.  The  governing  equation  on  the  unit  cube  fl  is 

jut-Au  =  +  20e-(“2(x-x2)^+(t-i2)")^  x  €  0,0  <  t, 

|  it(;c,0)  =  0,  X  fe  $2, 

with  homogeneous  Neumann  boundary  conditions  on  all  the  sides  except  the  bottom 
whore  a  homogeneous  Dirichlet  condition  is  imposed."  We  choose  ai  =  50,  tt2  =  10, 
ti  =  1,  1-2  =  10,  .7.1  =  (0.125,0.125,0.125)  and  .rj  -  (q.';j5, 0.5, 0.75).  The  quantity  of 
interest  is  the  time  average  of  the  solution  value. 

We  use  a  coarse  discretization  consisting  of  an  8  x  8  x  8  uniform  hexahedral  mesh 
and  time  step  of  0.1.  With  ATOL  =  0.00010044  and  Nmeix  =  50000  we  show  the 
block  information  for  the  memory  bound  and  correlation  strategies  respectively  in 
Table  4.2  and  Table  4.3.  The  algorithms  lead  to  significantly  different  block  meshes. 
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Block 

n i 

Tfc 

#  vertices 

■ff  hexahedra 

1 

0 

1.1 

59465 

54125 

2 

1.1 

1.2 

63112 

57772 

3 

1.2 

2.4 

45359 

40958 

4 

2.4 

11.9 

12383 

10165 

5 

11.9 

14.9 

2029 

1478 

Table  4.2 


Problem  (4-4)-  Blocks  resulting  from  the  memory  bound  strategy. 


Block 

n-i 

n 

if  vertices 

#  hexahedra 

1 

0 

1.1 

63112 

57772 

2 

1.1 

1.2 

63112 

57772 

3 

1.2 

1.6 

45359 

40958 

4 

1.6 

2.5 

9611 

8037 

5 

2.5 

2.9 

1968 

1436 

6 

2.9 

8.5 

966 

652 

7 

8.5 

9 

2617 

1926 

8 

9 

10.8 

12651 

10382 

9 

10.8 

11.3 

7363 

5860 

10 

11.3 

12.6 

3139 

2360 

11 

12.6 

14.9 

729 

512 

Table  4  3 

Problem  (4-4)-  Blocks  resulting  from  the  correlation  strategy. 


The  correlation  strategy  chooses  many  more  blocks,  but  many  of  the  blocks  have  very 
low  numbers  of  elements. 

We  show  planar  slices  near  xi  and  X2  of  the  meshes  for  Blocks  1  and  3  in  Fig.  4.5. 
For  comparison,  we  show  planar  slices  perpendicular  to  the  x-axis  near  xi  and  X2  of 
the  meshes  for  blocks  constructed  using  the  two  strategies  in  Fig.  4.6.  Both  strategies 
result  in  similar  meshes  near  X2  at  time  t  =  10.  However  at  t  =  8.8,  the  correlation 
strategy  leads  to  coarse  meshes  that  are  not  produced  by  the  memory  V)ound  strategy. 
The  mesh  resulting  from  the  memory  bound  strategy  retains  the  refinement  resulting 
from  the  earlier  perturbation  near  xi  at  t  =  1. 

4.5.  Example  Five:  Periodic  motion  in  a  convection-dominated  flow. 
This  example  has  a  heat  source  with  a  forced  oscillating  convective  term  within  the 
unit  cube  to  produce  an  “orbiting”  zone  of  perturbation.  The  governing  equation 
is 


{W(  +  ft  ■  Vw  -  Aji  =  /,  X  e  $2, 0  <  f  <  1, 
it(i,  t)  =  0,  X  e  dil,  0  <  f  <  1,  (4.5) 

u(x,0)  =  0,  X  €  n, 

with  >3  =  (20(cos(7rt)sin(27r<),sin(7rt)sin(27rt).cos(27rt))  and  /(x)  =  e-50(xj+xj  fxj) 
The  (juantity  of  interest  is  the  time  average  value.  The  coarse  discretization  used 
4913  vertices  and  at  time  step  of  0.01.  The  blocks  constructed  by  the  memory-bound 
strategy  using  ATOL  —  0.00044  and  Nniax=50000  arc  described  in  Table  4.4. 
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FlO.  4.5.  Problem  (4'4)-  Memory  bound  strategy.  Slice.s  through  the  mesh  perpendicular  to  the 
X-axis.  Upper  left:  Slice  near  xi  at  t  =  I  (block  I).  Upper  right:  Slice  near  xz  at  t  =  I  (block  1). 
Lower  left:  Slice  near  n  at  t  =  10  (block  4).  Lower  right:  Slice  near  xz  at  t  =  10  (block  4). 


Flo.  4.6.  Problem  (4.4) ■  Slices  through  the.  mesh  perpendicular  to  the  x- axis.  Left:  Correlation 
strategy.  Slice  near  xz  at  t  =  10  (block  8).  Middle:  Correlation  strategy.  Slice  near  xi  at  t  =  8.8 
(block  7).  Right:  Memory  bound  strategy.  Slice  near  ij  at  t  =  8.8  (block  4). 
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Block 

n 

T'h+l 

#  vertices 

#  hexahedra 

1 

0 

0.09 

58799 

51066 

2 

0.09 

0.15 

58424 

50289 

3 

0.15 

0.27 

58393 

50359 

4 

0.27 

0.61 

59102 

50744 

5 

0.61 

0.99 

28395 

23388 

TaU1,E  4.4 

Problein  (4-5).  Blocks  resulting  from  the  memory  bound  strategy. 


We  provide  “slices”  through  the  mesh  that  are  perpendicular  to  the  a;-axis  at 
X  =  0.5  for  four  representative  times  in  Fig.  4.7. 

5.  Conclusions.  In  this  paper,  we  consider  adaptive  algorithms  for  evolution 
problems  that  use  a  sequence  of  “blocks”  in  time  which  employ  fixed,  non-uniform 
space  meshes.  Blockwise  adaptive  algorithms  provide  a  way  to  balance  the  goal  of 
achieving  desired  accuracy  using  discretizations  with  relatively  few  degrees  of  freedom 
with  the  computational  overhead  associated  with  load  balancing,  re-meshing,  matrix 
reassembly  and  error  estimation.  Block  adaptive  algorithms  achieve  this  balance  by 
minimizing  the  number  of  mesh  changes.  However,  a  major  issue  is  determining  block 
discretizations  from  coarse  scale  solution  information  that  achieve  the  desired  accuracy 
and  efficiency.  We  describe  several  strategies  to  achieve  this  goal  using  adjoint-based 
a  posteriori  error  estimates.  Wo  demonstrate  the  behavior  of  the  proposed  algorithms 
as  well  as  several  technical  issues  in  a  set  of  examples. 
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Appendix  A.  Detailed  description  of  a  block  adaptive  algorithm. 


The  notation  used  in  our  block  adaptive  algorithm  is  as  follows. 


1.  Ntimestep  =  current  number  of  time  steps 

2.  Nelem(j)  -  number  of  .space  elements  in  timestep  j,  i.e.,  for  t.  € 

3.  Ntimestep.children(j)  =  number  of  subintervals  into  which  timestep  j  is 
to  be  divided 

4.  Nelem.children(i,  j)  --  number  of  subelemeuts  into  which  finite  element  i  is 
to  be  divided  in  timestep  j 

5.  The  6th  “block”  is  time  interval  [7b_i, Tb]  =  [tb.o, 4,ArJ 

6.  The  6th  “block”  comprises  timesteps  jb-i, ■  •  -Jb,  i-e->  Nb  -  jb  -  jb-\,  h,o  = 
ijb  1  ^iid  tb.jVi,  —  tj(,  • 

7.  block(i,b)  =  number  of  intervals  the  parent  element  i  will  be  divided  into 
on  block  6. 

8.  Nelem.block(b)  =  number  of  elements  in  block  b. 

9.  We  use  the  Matlab  colon  operator  :  to  denote  the  full  row  or  column. 

10.  The  parameter  0  governs  how  often  a  mesh  is  coarsened;  0  w  10  works  well. 
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Algorithm  A.l  A  memory-bound  strategy 

1:  Input  error  tolerance  TOL,  maximum  number  of  elements  in  any  block  Nmax,  the 

initial  coarse-scale  discretization  for  the  forward  problem,  and  the  coarse-scale 
discretization  for  the  adjoint  problems 
2:  Solve  forward  problem  (1.1)  for  U  on  [0,  T] 

3:  Project  forward  solution  onto  coarse-scale  adjoint  problem  mesh 
4:  Solve  adjoint  problem  (2.4)  on  coarse  scale  mesh  and  compute  E{U) 

5:  Compute  LATOL,  JEt  (3. 6), (3. 5) 

6:  for  j  =  1, . . . ,  Ntimesteps  do 
7:  Compute  Ntimestep_children(j )  (3.13) 

8:  for  f  =  1, . . . ,  Nelem(  j )  do 

9:  Compute  Nelem.childrenCi.j)  (3.16) 

10:  end  for 

11:  end  for 

12:  Ntimesteps  <—  Ntimestep_cliildren(j) 

13:  Each  subinterval  of  inherits  Nelem_children(i,  j) 

M:  6=1,  To  =  0,  1\  =  ki,  jo  r=  1,  J  =  2 

15:  block(;,b)  <— Nelem_children(;,  1) 

16:  Nelem_block(b)  <—  ^jblock(i,b) 

17:  while  Tfc  <  T  do 

18:  while  Nelem,.block(b)  <  Nmax  and 

19:  Nelem-block(b)  <0  x  Nelem-children(i,  j)  do 

20:  jb  i-  j 

21:  Tb  Tb  b  kj 

22:  block(:,b)  <- max[block(;,b),Nelem.children(:,  j)] 

23:  Nelem.block(b)  =  block(i,  b) 

24:  j  <-  j  +  1 

25:  end  while 

26:  6  e-  6  -4-  1 

27:  end  while 

28:  for  1  =  1, ...  ,6  do 

29:  Compute  new  mesh  for  block  b 

30:  Optional:  Estimate  projection  error  and  correct  predicted  meshes  if  necessary 

31:  end  for 

32:  for  i  =  1, . . . ,  6  do 

33:  Solve  forward  problem  on  block  b  for  U 

34:  Project  U  onto  mesh  for  block  6-1-1 

35:  Optional:  Compute  projection  error  between  blocks  and  correct  meshes 

36:  end  for 


To  implement  the  correlation-based  strategy,  we  alter  the  block  selection  criteria 
(X!block(6)  <  Nmax)  with  a  step  which  accepts  a  block  if  block(:,6)  is  correlated  to 
Nelem_children(;,  j)  and  Nelem.block(b)  is  less  than  Nmax. 

The  algorithm  assumes  that  the  blocks  are  always  generated  (even  on  repeat  solve 
c:yc:lcs)  using  the  coarse  mesh  as  a  base.  The  algorithm  may  lie  (;asily  modified  to 
work  recursively  on  the  blocks.  It  may  also  be  modified,  with  a  little  more  care,  to 
allow  merging  and  splitting  of  blocks  during  repeated  solves. 
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